1. Financial Data & Returns
1.1 Limit Order Book
Go to https://lobsterdata.com/info/DataSamples.php and download the level 10 data for Amazon stock (AMZN). Use the provided sample R code to process the data; you can find a description of the resulting R data-frames (i.e. message & order book file) generated by the sample R code here. Use the data to create a daily plot of the evolution of the limit order book at 5-minute intervals
# load lobster demo data in R
source("./res/ps1/LOBSTER_demo.R", chdir = TRUE)
## [1] "No trading halts detected."
# trading happens from 9:30pm to 4:00pm
# in term of seconds, this is from 34,200 to 57,600
# (note that the time variable in the first column of data is the time of day in seconds)
times = seq( from=34200, to=57600, by = 5*60 ) # 5 minute intervals
timeindex=rep(0,length(times)) # index to hold row numbers of 5min intervals
ASKP=ASKV=BIDP=BIDV=matrix(0,length(times),10) # price and volume at 5min intervals
for(i in 1: length(times)){
timeindex[i] = which.min( abs(dataM$Time-times[i]) )
ASKP[i,] = as.numeric(dataOB[timeindex[i],(0:9)*4+1])
ASKV[i,] = as.numeric(dataOB[timeindex[i],(0:9)*4+2])
BIDP[i,] = as.numeric(dataOB[timeindex[i],(0:9)*4+3])
BIDV[i,] = as.numeric(dataOB[timeindex[i],(0:9)*4+4])
}
MIDP=(ASKP[,1] + BIDP[,1])/2 # midpoint of best bid & ask
ylim=range( cbind(ASKP,BIDP) ) + c(-1,+1)*10000
xlim=range( times )
plot(-1,-1, xlim = xlim, ylim = ylim, xlab="time of day (sec)", ylab='Price ($1/10,000)')
# plots of 10 best bid/ask prices at 5min intervals
# volume at each price level is represented by size of points
for(i in 1:length(times)){
xloc = times[i] + 10*(1:10)
points( xloc, ASKP[i,], col=2, pch='.', cex=(ASKV[i,])^(1/3) )
points( xloc, BIDP[i,], col=3, pch='.', cex=(BIDV[i,])^(1/3) )
}
lines(times, MIDP)
# grid
abline(h=seq(2200000, 2300000, by=10000), lty=3, col='lightgrey')
abline(v=(10:16)*60*60, lty=3, col="lightgrey")
legend('topright', title="Volume", c("100","250",'1000'), pch='.', pt.cex=c(100,250,1000)^(1/3))
1.2 Daily Return
Download the closing and adjusted closing daily prices for Apple Inc.
(symbol AAPL) for the previous year (Jan 1, 2022 to Dec 31, 2022). Use
the function get.hist.quote() from the tseries package.
Plot the closing and adjusted-closing price series on the same plot. Do you see any differences?
library(tseries)## Registered S3 method overwritten by 'quantmod': ## method from ## as.zoo.data.frame zoo# download prices S.cl=get.hist.quote("AAPL", start='2022-01-01', end='2022-12-31', quote='Close')## time series starts 2022-01-03 ## time series ends 2022-12-30S.adj=get.hist.quote("AAPL", start='2022-01-01', end='2022-12-31', quote='AdjClose')## time series starts 2022-01-03 ## time series ends 2022-12-30plot( S.cl ) lines( S.adj, col = 2 ) legend("topright", lwd=2, col = 1:2, c("Close", "AdjClose"))Calculate and plot the net returns based on the adjusted closing price.
R = diff(S.adj) / timeSeries::lag(S.adj,-1) plot(R)# # Alternatively, you can use the function quantmod::Delt() # library(quantmod) # R = Delt( S.adj ) # # Note: for log-returns you can use either of: # r = Delt( S.adj, type = "log" ) # r = diff( log( S.adj ) ) # # Also note that: R = exp( r ) - 1 <=> r = log( R + 1 )Apple paid a dividend of $0.23 on Nov 4, 2022. Calculate the net return of Nov 4, using both the adjusted closing prices (designated way), and the actual closing prices, where you manually adjust for the dividend (i.e., you add the dividend amount to the Nov 4 close price before calculating the log-return). Are the two returns equal?
The return calculated based on the AdjClose prices is:
R[as.Date("2022-11-04")]## 2022-11-04 ## -0.001947376Note:
get.hist.quote()creates time series that are irregularly-spaced (because of weekday trading days, holidays, etc) and cannot be represented by the base-R classts. Therefore, the function outputs an object of classzoo, which is a precursor and special case of classxts. For more details on these classes and on working with irregularly spaced time series in R, see Manipulating Time Series Data in R with xts & zooThe return calculated based on the dividend correction is:
Div = .23 Pt = S.cl[as.Date("2022-11-04")] Pt_1 = as.numeric( S.cl[as.Date("2022-11-03")] ) # convert to numeric so that time indexes don't clash ( ( Pt + Div) - Pt_1 ) / ( Pt_1 )## 2022-11-04 ## -0.001944124This is slightly different than the return based on the adjusted closing price. The reason for that is how Yahoo! Finance implements its dividend adjustment, which is multiplicative rather than additive, in order to make it apply cumulatively to all previous prices (hence the difference in Close and AdjClose values before Nov 4 2022). For more details on how Yahoo! Finance calculates adjusted close prices see here.
Consider a simple momentum strategy whereby:
- you buy APPLE stock for 1 day, if its price increased by more than 2% on the previous day, and
- you (short) sell APPLE stock for 1 day, if its price decreased by more than 2% on the previous day. Calculate and plot the daily and cumulative net returns of this strategy.
mom.ret = sign( timeSeries::lag(R, -1) ) * R # ts of momentum strategy daily returns cum.mom.ret = cumprod( 1 + mom.ret ) - 1 # cumulative returns ( have to multiply/compound net returns ) plot(cbind(mom.ret,cum.mom.ret), main="", ylab=c("Daily Returns","Cumulative Returns") , xlab='Time')Find the maximum
drawdownof the strategy, i.e. the biggest drop in cumulative net returns from their highest point. You can use themaxdrawdown()function from thetseriespackage.MDD = maxdrawdown(cum.mom.ret) MDD## $maxdrawdown ## [1] 0.5041224 ## ## $from ## [1] 61 ## ## $to ## [1] 230plot(cum.mom.ret,ylab="Cumulative Returns" , xlab='Time') lines( cum.mom.ret[c(MDD$from[1],MDD$to[1])], col=2 ) abline(h=0,lty=2)
1.4 Probability calculations
Assume that the daily log returns on a stock with current price \(S_0 = 100\) are i.i.d. with a Normal distribution with mean 0.0003 and standard deviation 0.01.
Calculate the probability that the price will be above $120 after 50 days.
We have \[ \begin{aligned} P( S_{50} > 120 ) &= 1 - P( S_0\exp\{r_{1-50}\} \le 120 ) \\ &= 1 - P( r_{1-50} \le \ln(120/100) )\\ &= 1 - P( Z \le \frac{ \ln(1.2) - .0003\times 50 }{ 0.01\times \sqrt{50} } ) \end{aligned} \] which gives
1 - pnorm( (log(1.2) - .0003*50)/(.01*sqrt(50)) )## [1] 0.008983824Calculate the probability that the price will increase (i.e., have +ve return) in 30 or more of the next 50 days. Note: You will have to use (pnorm(), pbinom()) to find the numerical values of these probabilities
The probability the price will increase on any given day is \[ \begin{aligned} P(S_t > S_{t-1}) &= P(S_t/S_{t-1}>1) \\ & = P( \ln(S_t/S_{t-1}) > \ln(1) ) \\ & = P( r_t > 0 ) = P( Z > \frac{0 - 0.0003}{0.01} ) \\ & = \Phi(0.03) \end{aligned} \] which is
p = pnorm(0.03)The probability that the price will increase in at least 30 out of 50 days is equal to the probability of a Binomial(\(n=50,p=P(S_t > S_{t-1})\)) RV \(X\) taking a value greater or equal to 30, i.e. \(P(X\ge 30 ) = 1 - P(X < 30) = 1-P(X \le 29)\), which is
1 - pbinom( 29, 50, p)## [1] 0.1346917or, equivalently,
pbinom( 29, 50, p, lower.tail = F)## [1] 0.1346917
2. Univariate Return Modelling
2.4 Return Distributions
Download from Yahoo! Finance the daily adjusted closing prices of Loblaw Companies Limited (L.TO), from Jan-1-2015 to Dec-31-2022:
S = get.hist.quote( "L.TO", start='2015-01-01', end='2022-12-31', quote='AdjClose', drop = TRUE)
## time series starts 2015-01-02
## time series ends 2022-12-30
Calculate the daily net returns, and plot their time-series plot, their sample auto correlation plot, and their Normal QQ plot. Which of the return stylized facts can you observe?
R = quantmod::Delt(S)[-1] # Remove first NA return plot(R)acf( as.numeric(R) )qqnorm(R, datax = TRUE); qqline(R, datax = TRUE)The returns exhibits volatility clustering (time series plot), are uncorrelated (ACF plot), and have much heavier tails than Normal (Normal QQ-plot).
Fit a univariate t-distribution to the returns and report the resulting values (use the fitdistr function from the MASS library). Find the mean and variance of the fitted t-distribution, and compare them to the sample mean and variance of the returns.
library(MASS) t_est = fitdistr(R, "t")$estimate # estimated parameters are: (mu, sigma, v) # Normal-chi mixture model: N(mu, sigma) * sqrt( v / chi2(v) ) mu = t_est[1] sigma = t_est[2] v = t_est[3] # t mean & variance mean.t = mu var.t = sigma^2 * ( v/(v-2) ) # the sample moments are mean.norm = mean(R) var.norm = var(R)The variance of the fitted t distribution (1.4635153^{-4}) is close to the sample variance (1.3827259^{-4}), while the mean of the t distribution (3.2241245^{-4}) is lower than the sample mean (5.661532^{-4}). Differences can occur in the parameters because the t distribution can absorb extreme values.
We will now compare different return distribution approaches in an practical investment setting. Assume you invest all of your wealth in L.TO for 4 years (4 × 252 = 1008 days). Simulate 5,000 iterations of 1008 daily returns from the following models:
- Returns are i.i.d. Normal with parameters equal to the sample mean and variance (i.e., the Normal MLE estimates). in 1.b
- Returns are i.i.d. t with the parameters identified in the previous part.
- (Optional) Returns follow a GARCH(1,1) model: use the garchFit() function to fit the model and the garchSim function to simulate from it (both in library fgarch).
Calculate the final wealth and the maximum drawdown of each iteration, and create a histograms of the two metrics from each model (2 × 3 = 6 histograms in total). What do you observe?
n = 5000; m = 1008 set.seed(1234567890) library(fGarch)## NOTE: Packages 'fBasics', 'timeDate', and 'timeSeries' are no longer ## attached to the search() path when 'fGarch' is attached. ## ## If needed attach them yourself in your R script by e.g., ## require("timeSeries")### i # create n x m matrix of simulated paths (where each row = a path) R.norm = matrix( rnorm( n*m, mean.norm, sqrt(var.norm) ), nrow = n, ncol = m ) # apply cumulative product to each row, to get the cumulative net return path Rcm.norm = t( apply( 1+R.norm, MARGIN = 1, FUN = cumprod) - 1 ) # MARGIN=1 for rows, 2 for cols # final returns are the last column/value in each path Rfn.norm = Rcm.norm[,m] # calculate maximum drawdown for each path mdd = function(R){ return(tseries::maxdrawdown(R)$maxdrawdown) } MDD.norm = apply(Rcm.norm, MARGIN = 1, FUN = mdd ) ### ii # similarly for t distribution R.t = matrix( rnorm( n*m, mu, sigma ) * sqrt( v / rchisq( n*m , df = v) ), nrow = n, ncol = m ) Rcm.t = t( apply( 1+R.t, MARGIN = 1, FUN = cumprod) - 1 ) Rfn.t = Rcm.t[,m] MDD.t = apply(Rcm.t, MARGIN = 1, FUN = mdd ) ### iii GARCH.fit=garchFit(~garch(1,1), data=R, trace = FALSE) # fit GARCH(1,1) model GARCH.param=GARCH.fit@fit$coef # fitted coefficients # define GARCH model specification for simulation GARCH.spec=garchSpec(model=list(mu=GARCH.param['mu'], omega=GARCH.param['omega'], alpha=GARCH.param['alpha1'], beta=GARCH.param['beta1'] )) R.garch = matrix( 0, nrow = n, ncol = m ) for(i in 1:n){ R.garch[i,]=as.numeric(garchSim(GARCH.spec,m)) } Rcm.garch = t( apply( 1+R.garch, MARGIN = 1, FUN = cumprod) - 1 ) Rfn.garch = Rcm.garch[,m] MDD.garch = apply(Rcm.garch, MARGIN = 1, FUN = mdd )par(mfrow = c(2,3)) hist(Rfn.norm) hist(Rfn.t) hist(Rfn.garch) hist(MDD.norm) hist(MDD.t) hist(MDD.garch)The final returns for t-distr. and GARCH gave lower negative values, but are comparable to those of Normal. In terms of maximum drawdown, the GARCH values have a smaller range, but are otherwise also comparable. Not that for for both metrics we combine/compound many returns, so the effects of heavy tails is suppressed.
2.5 Heavy Tail Distribution Results & EVT
In this question you will verify empirically (using simulation) some basic results from heavy tailed distributions and extreme value theory.
Generate \(n=5000\) i.i.d. Cauchy variates, and calculate the cumulative sample average \(\bar{X}_k=\frac{1}{k} \sum_{i=1}^k X_i\) versus \(k\), for \(k=1, \ldots, n\). If the WLLN held for this distribution, you would expect to see the average converge to 0 . Does the plot behave like that? (Hint: the Cauchy distribution is the same as a \(t\) with 1 degree of freedom, use this fact to generate Cauchy variates using the function
rt().)set.seed(1234567890) n=5000 plot( cumsum( rt(n , df = 1) ) / (1:n), type = "l"); abline(0,0, lty =2)The sample average of the Cauchy distribution does not converge to a its mean (0), as there will always be large shifts caused by extreme values. You can compare this with the plot of the average of i.i.d. Normals, for which the WLLN applies, and which shows the expected convergence.
plot( cumsum( rnorm(n) )/(1:n), type = "l"); abline(0,0, lty =2)Generate \(n=5000\) pairs \(\left(X_i, Y_i\right)\) of independent Cauchy variates, and calculate their mean \(W_i=\left(X_i+Y_i\right) / 2\). The mean should also be Cauchy distributed, according to Q2c (i.e., Cauchy is a stable distribution). Create a QQ-plot comparing the sample quantiles from your simulation to the theoretical quantiles from a Chauchy distribution; does the plot confirm the theoretical result. (Hint: use the function
qt(ppoints(n, df=1)to find the theoretical quantiles of the Cauchy distribution.)X = rt(n, df=1) Y = rt(n, df=1) W = (X + Y)/2 x = sort(W) #sample quantiles y = qt( ppoints(n), df = 1) # theoretical quantiles plot( x, y , xlab= "Sample Quantiles", ylab="Theoretical Quantiles"); abline(0,1)The QQ plot shows good agreement of the Cauchy average quantiles vs the theoretical quantiles of the Cauchy distribution, as expected by the theoretical result. (The deviations at the very ends of the tails, where there are fewer values, are expected to some extent for extreme heavy-tail distributions like the Cauchy).
The 1st Extreme Value Theorem (Fisher-Tippett-Gnedenko) states that the rescaled maximum \(\frac{M_n-b_n}{a_n}\) of n i.i.d. RVs can only converge to one of three types of distributions. Repeat the following simulation experiment:
Generate \(m=1000\) values of the maximum of \(n=100\) i.i.d. RVs, i.e., generate \(M_n(j)=\max \left\{X_i(j), \ldots, X_n(j)\right\}\) for each iteration \(j=1, \ldots, m\).
Rescale the simulated maxima according to appropriate \(a_n, b_n\), i.e., calculate \(Y_j=\frac{M_n(j)-b_n}{a_n}\) for \(j=1, \ldots, m\).
Plot the histogram of the rescaled maxima \(Y_j\), overlaid with the density of the corresponding theoretical limiting distribution (Weibull, Gumbel, or Fréchet). for each of the following three distributions and associated scaling parameters:
\(X_i \sim\) i.i.d. Uniform(0, 1) with \(b_n=1\) and \(a_n=1 / n\); which converges to Weibull \((\) with \(\alpha=1\) ).
\(X_i \sim\) i.i.d. Normal(0,1) with \(b_n=\Phi^{-1}(1-1 / n)\) and \(a_n=\frac{1}{n \phi\left(b_n\right)}\); which converges to Gumbel. (Where \(\Phi / \phi\) are the std Normal CDF/PDF.)
\(X_i \sim\) i.i.d. Cauchy with \(b_n=0\) and \(a_n=n / \pi\); which converges to Fréchet (with \(\alpha=1\) ).
Does the theorem seem to hold?
n=100; m=1000 ### i # Generate iid RVs in (m x n) matrix U = matrix( runif(m*n), m, n) # Calculate normalized maxima over rows b_n = 1 a_n = 1/n Mn = (apply(U, 1, max)-b_n)/a_n # Plot histogram hist(abs(Mn), probability = T, main = "i.") # overlay theoretical distribution x=seq(0,10,.01) lines(x, dexp(x), col = 2) legend( "topright", lwd = 2, col = 1:2, c("Sample hist.","Weibull"))### ii Z = matrix( rnorm(m*n), m, n ) b_n = qnorm( 1 - 1/n ) a_n = 1 / (n * dnorm(b_n) ) Mn = (apply(Z, 1, max) - b_n ) / (a_n) hist(Mn, probability = T, main = "ii.") x=seq(-5,5,.01) lines(x, exp( - x - exp(-x)), col = 2) legend( "topright", lwd = 2, col = 1:2, c("Sample hist.","Gumbel"))### iii X= matrix( rt(m*n,df=1), m, n ) Mn = (apply(X, 1, max) ) / (n/pi) hist(Mn, probability = T, main = "iii." ) x=seq(0,150,.01) lines(x, 1/x^2 * exp(- 1/x), col = 2) legend( "topright", lwd = 2, col = 1:2, c("Sample hist.","Frechet"))In cases i. and ii. the convergence to the theoretical distribution is evident, but in iii. the histogram is not as informative because of the extreme values involved. For the last case, one can plot the histogram & density of the log-transformed normalized maxima, which actually show the convergence. (Note: the log-transformation of the Frechet(\(\alpha=1\)) is the Gumbel distribution.)
hist( log(Mn), probability = T, main = "iii. log-transform" ) x=seq(-2,15,.01) lines(x, exp( -x-exp(-x)), col = 2) legend( "topright", lwd = 2, col = 1:2, c("log-Sample hist.","log-Frechet"))
3. Multivariate Return Modelling
3.4 Equity Returns
This section uses the data set berndtInvest on the
book’s web site and taken originally from R’s fEcofin
package. This data set contains monthly returns from January 1, 1987, to
December 1, 1987, on 16 equities. There are 18 columns. The first column
is the date and the last is the risk-free rate.
In the lab we will only use the first four equities. The following code computes the sample covariance and correlation matrices for these returns.
berndtInvest = read.csv("./res/datasets/berndtInvest.csv")
Berndt = as.matrix(berndtInvest[, 2:5])
cov(Berndt)
## CITCRP CONED CONTIL DATGEN
## CITCRP 0.006556442 0.0010933247 0.006052187 0.0055091027
## CONED 0.001093325 0.0025272604 0.000795463 0.0006124665
## CONTIL 0.006052187 0.0007954630 0.022710259 0.0050502504
## DATGEN 0.005509103 0.0006124665 0.005050250 0.0162664367
cor(Berndt)
## CITCRP CONED CONTIL DATGEN
## CITCRP 1.0000000 0.2685901 0.4959835 0.5334584
## CONED 0.2685901 1.0000000 0.1049987 0.0955237
## CONTIL 0.4959835 0.1049987 1.0000000 0.2627578
## DATGEN 0.5334584 0.0955237 0.2627578 1.0000000
If you wish, you can also plot a scatterplot matrix with the following R code.
pairs(Berndt)
Q1
Suppose the four variables being used are denoted by \(X_1, \ldots, X_4\). Use the sample covariance matrix to estimate the variance of \(0.5 X_1+0.3 X_2+\) \(0.2 X_3\).
COV = var(Berndt) # variance-covariance matrix of returns
w = c(.5, .3, .2, 0) # vector of weights
t(w) %*% COV %*% w # variance of linear combination of returns (portfolio)
## [,1]
## [1,] 0.004408865
Fit a multivariate- \(t\) model to
the data using the function cov.trob in the MASS package.
This function computes the MLE of the mean and covariance matrix with a
fixed value of \(\nu\). To find the MLE
of \(\nu\), the following code computes
the profile log-likelihood for \(\nu\).
library(MASS) # needed for cov.trob
library(mnormt) # needed for dmt
df = seq(2.5, 8, 0.01)
n = length(df)
loglik_profile = rep(0, n)
for (i in 1:n) {
fit = cov.trob(Berndt, nu = df[i])
mu = as.vector(fit$center)
sigma = matrix(fit$cov, nrow = 4)
loglik_profile[i] = sum(log(dmt(Berndt, mean = fit$center, # mt for multivar t
S = fit$cov, df = df[i])))
}
plot(loglik_profile, type = 'l')
Q2
Using the results produced by the code above, find the MLE of \(\nu\) and a 90% profile likelihood confidence interval for \(\nu\). Include your R code with your work. Also, plot the profile log-likelihood and indicate the MLE and the confidence interval on the plot.
Section 7.13.3 demonstrates how the MLE for a multivariate t-model
can be fit directly with the optim function, rather than
profile likelihood.
plot(df, 2 * loglik_profile, type = "l", cex.axis = 1.5, # plot 2*log(likelihood)
cex.lab = 1.5, ylab = "2 * log(likelihood)", lwd = 2)
abline(h = 2 * max(loglik_profile) )
# plot maximum
max_ind = which.max(loglik_profile) # Find index of maximum log(likelihood)
points( df[max_ind], 2*loglik_profile[max_ind], pch=15, cex=2)
text( df[max_ind], 2*loglik_profile[max_ind], pos=1, labels = "MLE")
# plot 90% CI level
abline(h = 2 * max(loglik_profile) - qchisq(0.9, 1))
# find indexes of lower & upper bounds of 90% CI
lower_ind = head( which( loglik_profile > max(loglik_profile) - qchisq(0.9, 1 )/2 ), 1 )
upper_ind = tail( which( loglik_profile > max(loglik_profile) - qchisq(0.9, 1 )/2 ), 1 )
abline(v = df[lower_ind] )
abline(v = df[upper_ind] )
points( df[lower_ind], 2*loglik_profile[lower_ind], pch=16, cex=2) # plot lower bound
text( df[lower_ind], 2*loglik_profile[lower_ind], pos=1, labels = "lower 90% \n CI bound")
points( df[upper_ind], 2*loglik_profile[upper_ind], pch=16, cex=2) # plot upper bound
text( df[upper_ind], 2*loglik_profile[upper_ind], pos=1, labels = "upper 90% \n CI bound")
print( rbind( c( "lower", "max", "upper"), df[ c(lower_ind, max_ind, upper_ind)] ))
## [,1] [,2] [,3]
## [1,] "lower" "max" "upper"
## [2,] "3.23" "4.55" "6.59"
Q3
The following code generates and plots four bivariate samples. Each sample has univariate marginals that are standard t3-distributions. However, the dependencies are different.
library(MASS) # need for mvrnorm
par(mfrow=c(2,2))
N = 2500
nu = 3
set.seed(5640)
cov=matrix(c(1, 0.8, 0.8, 1), nrow = 2) # rho = 0.8 (correlated)
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w, w) # scale by same factor (tail dep)
plot(x, main = "(a)")
set.seed(5640)
cov=matrix(c(1, 0.8, 0.8, 1),nrow = 2) # rho = 0.8
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w1 = sqrt(nu / rchisq(N, df = nu))
w2 = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w1, w2) # scale by dif factors
plot(x, main = "(b)")
set.seed(5640)
cov=matrix(c(1, 0, 0, 1), nrow = 2) # rho = 0
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w1 = sqrt(nu / rchisq(N, df = nu))
w2 = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w1, w2) # scale by dif factors
plot(x, main = "(c)")
set.seed(5640)
cov=matrix(c(1, 0, 0, 1), nrow = 2) # rho = 0
x= mvrnorm(N, mu = c(0, 0), Sigma = cov)
w = sqrt(nu / rchisq(N, df = nu))
x = x * cbind(w, w) # scale by same factor
plot(x, main = "(d)")
Which sample has independent variates? Explain your answer.
- Samples (a) and (b) are (positively) correlated, because they both arise from (randomly scaled) correlated Normals (with \(\rho=0.8\)).
- Sample (d) exhibits tail-dependence, because the uncorrelated/independent Normals are scaled by the same random (\(\chi^2\)-based) factor.
- Only sample (c) is independent, because the two independent Normal variates are scaled by independent random factors.
Q4
Which sample has variates that are correlated but do not have tail dependence? Explain your answer.
- Sample (b) has correlated variates (because it is based on correlated Normals), but does not have tail-dependence, because its variates were scaled by independent random factors. The plot shows that extreme values seem to occur independently in the two dimensions (i.e., extreme values are either horizontally or vertically extreme, but not both.)
Q5
Which sample has variates that are uncorrelated but with tail dependence? Explain your answer.
- Sample (d) has uncorrelated variates (based on uncorrelated Normals) but with tail-dependence, because they are scaled by the same random factor. The plot shows that extreme values tend to occur simultaneously in both dimensions.
3.5 Copulas
Run the R code that appears below to generate data from a copula.
library(copula)
# define copula obj
cop_t_dim3 = tCopula(dim = 3, param = c(-0.6,0.75,0), dispstr = "un", df = 1)
set.seed (5640)
# generate random sample from copula obj above
rand_t_cop = rCopula(n = 500, copula = cop_t_dim3)
pairs(rand_t_cop)
cor(rand_t_cop)
## [,1] [,2] [,3]
## [1,] 1.0000000 -0.54999514 0.70707296
## [2,] -0.5499951 1.00000000 -0.06538499
## [3,] 0.7070730 -0.06538499 1.00000000
dispstr is dispersion structure. It is the type of the
symmetric positive definite matrix characterizing the elliptical copula.
Currently available structures are “ex” for exchangeable, “ar1” for
AR(1)AR(1), “toep” for Toeplitz (toeplitz), and “un” for
unstructured.
Q1
Consider the R code above.
What type of copula model has been sampled? Give the copula family, the correlation matrix, and any other parameters that specify the copula.
It is a t copula. The correlation matrix is “unstructured” meaning that it is an arbitrary correlation matrix. In fact, the correlation matrix is
\[ \left(\begin{array}{ccc} 1 & -0.6 & 0.75 \\ -0.6 & 1 & 0 \\ 0.75 & 0 & 1 \end{array}\right) \] The degrees of freedom parameter is 1.
Q2
Examine the scatterplot matrix (generated by line 6) and answer the questions below. Include the scatterplot matrix with your answer.
Components 2 and 3 are uncorrelated. Do they appear independent? Why or why not?
Components 2 and 3 would be uniformly scattered over the unit square if they were independent.
Clearly, the scatter is not uniform, so they do not appear independent.
Do you see signs of tail dependence? If so, where?
The non-uniformity mentioned in (a) is that there are more data in the corners, which shows that extreme values tend to occur together, although because of the zero correlation, a positive extreme value of one component is equally likely to be paired with a positive or negative extreme value of the other component.
What are the effects of dependence upon the plots?
The effects of tail dependence is the tendency of extreme values to pair. The negative correlation of components 1 and 2 shows in the concentration of the data along the diagonal from upper left to lower right. Positive extreme values in one component tend to pair with negative extreme values of ther other component.
The positive correlation of components 2 and 3 shows in the concentration of the data along the diagonal from lower left to upper right. Positive extreme values in one component tend to pair with positive extreme values of ther other component.
The nonzero correlations in the copula do not have the same values as the corresponding sample correlations. Do you think this is just due to random variation or is something else going on? If there is another cause besides random variation, what might that be? To help answer this question, you can get confidence intervals for the Pearson correlation: For example,
cor.test(rand_t_cop[,1],rand_t_cop[,3])will give a confidence interval (95 percent by default) for the correlation (Pearson by default) between components 1 and 3. Does this confidence interval include 0.75?The question is asking why cor(ran_t_cop) is different from the correlation params passed into the copula
cop_t_dim3.cor.test(rand_t_cop[,1],rand_t_cop[,3])## ## Pearson's product-moment correlation ## ## data: rand_t_cop[, 1] and rand_t_cop[, 3] ## t = 22.314, df = 498, p-value < 2.2e-16 ## alternative hypothesis: true correlation is not equal to 0 ## 95 percent confidence interval: ## 0.6603249 0.7483624 ## sample estimates: ## cor ## 0.707073The confidence interval is (0.6603, 0.7484) which does not include 0.75, which was passed into the copula object.
0.75 is the correlation between the t-distributed random variables that define the copula and need not be the same as the uniformly-distributed variables in the copula itself.
Q3
# create Gaussian copula obj
cop_normal_dim3 = normalCopula(dim = 3, param = c(-0.6,0.75,0),
dispstr = "un")
# define dist using mvdc (multivar dist constructed from copulas)
mvdc_normal = mvdc(copula = cop_normal_dim3, margins = rep("exp", 3),
paramMargins = list(list(rate=2), list(rate=3), list(rate=4)))
# generate random sample
set.seed(5640)
rand_mvdc = rMvdc(n = 1000, mvdc = mvdc_normal)
pairs(rand_mvdc)
par(mfrow = c (2,2))
for(i in 1:3) plot(density(rand_mvdc[,i]))
What are the marginal distributions of the three components in rand mvdc? What are their expected values?
The marginal distributions for all 3 components are Exponential with rate \(\lambda= 3\), which means that their mean value is \(1/\lambda = 1/3\).
Are the second and third components independent? Why or why not?
The 2nd and 3rd components are (marginally) independent, since they jointly (i.e., var2 and var3 combined, without var1) follow 2D Normal/Gaussian, and their correlation is 0.
Their correlation is given by the
copula::p2Pfunction, which specifies the correlation matrix based on the vector of distinct correlations; for this problem we had:# p2P() creates a matrix from a given vector of parameters. # P2p() creates a numeric vector from a given matrix p2P( c( -0.6, +0.75, 0 ) )## [,1] [,2] [,3] ## [1,] 1.00 -0.6 0.75 ## [2,] -0.60 1.0 0.00 ## [3,] 0.75 0.0 1.00which means that
- var1 & var2 are negatively correlated,
- var1 & var3 are positively correlated, and
- var2 & var3 are un-correlated
Note that these relations are reflected in the Exponential random variates that were generated. Also note that the conditional distribution var2 and var3 given var1 would NOT be independent (that’s why we have to specify that the “marginal” distribution of var2 & var3 is independent.)
4. Portfolio Theory
4.2
This section uses daily stock prices in the data set
Stock_Bond.csv that is posted on the book’s website and in
which any variable whose name ends with “AC” is an adjusted closing
price. As the name suggests, these prices have been adjusted for
dividends and stock splits, so that returns can be calculated without
further adjustments. Run the following code which will read the data,
compute the returns for six stocks, create a scatterplot matrix of these
returns, and compute the mean vector, covariance matrix, and vector of
standard deviations of the returns. Note that returns will be
percentages.
dat = read.csv("./res/datasets/Stock_Bond.csv", header = T)
prices = cbind(dat$GM_AC, dat$F_AC, dat$CAT_AC, dat$UTX_AC,
dat$MRK_AC, dat$IBM_AC)
n = dim(prices)[1]
returns = 100 * (prices[2:n, ] / prices[1: (n-1), ] - 1)
pairs(returns)
mean_vect = colMeans(returns)
cov_mat = cov(returns)
sd_vect = sqrt(diag (cov_mat))
Q1
Write an R program to find the efficient frontier, the tangency portfolio, and the minimum variance portfolio, and plot on “reward-risk space” the location of each of the six stocks, the efficient frontier, the tangency portfolio, and the line of efficient portfolios.
Use the constraints that \(−0.1 \leq w_j \leq 0.5\) for each stock. The first constraint limits short sales but does not rule them out completely. The second constraint prohibits more than 50% of the investment in any single stock.
Assume that the annual risk-free rate is 3% and convert this to a daily rate by dividing by 365, since interest is earned on trading as well as nontrading days.
M = length(mean_vect)
library(quadprog)
Amat = cbind(rep(1,M),mean_vect,diag(1,nrow=M),-diag(1,nrow=M))
muP = seq(.05,0.08,length=300)
sdP = muP
weights = matrix(0,nrow=300,ncol=M)
# quadratic programming problem
for (i in 1:length(muP)){
result = solve.QP(Dmat=cov_mat, dvec=rep(0,M), Amat=Amat,
c(1,muP[i],rep(-.1,M),rep(-.5,M)), meq=2)
sdP[i] = sqrt(2*result$value)
weights[i,] = result$solution
}
plot(sdP,muP,type="l",xlim=c(0,2.5),ylim=c(0,.1))
mufree = 3/365
points(0,mufree,cex=3,col="blue",pch="*")
sharpe = (muP-mufree)/sdP
ind = (sharpe == max(sharpe)) # locates the tangency portfolio
weights[ind,] # weights of the tangency portfolio
## [1] -0.091181044 -0.002910879 0.335318542 0.383714329 0.319484849
## [6] 0.055574204
lines(c(0,sdP[ind]),c(mufree,muP[ind]),col="red",lwd=3)
points(sdP[ind],muP[ind],col="blue",cex=3,pch="*")
ind2 = (sdP == min(sdP))
points(sdP[ind2],muP[ind2],col="green",cex=3,pch="*")
ind3 = (muP > muP[ind2])
lines(sdP[ind3],muP[ind3],type="l",xlim=c(0,.25),
ylim=c(0,.3),col="cyan",lwd=3)
text(sd_vect[1],mean_vect[1],"GM")
text(sd_vect[2],mean_vect[2],"F")
text(sd_vect[3],mean_vect[3],"UTX")
text(sd_vect[4],mean_vect[4],"CAT")
text(sd_vect[5],mean_vect[5],"MRK")
text(sd_vect[6],mean_vect[6],"IBM")
legend("topleft",c("efficient frontier","efficient portfolios",
"tangency portfolio","min var portfolio"),
lty=c(1,1,NA,NA),
lwd=c(3,3,1,1),
pch=c("","","*","*"),
col=c("cyan","red","blue","green"),
pt.cex=c(1,1,3,3)
)
Q2
If an investor wants an efficient portfolio with an expected daily return of 0.07%, how should the investor allocate his or her capital to the six stocks and to the risk-free asset? Assume that the investor wishes to use the tangency portfolio computed with the constraints \(−0.1 \leq w_j \leq 0\)., not the unconstrained tangency portfolio.
Let \(\omega\) be the weight for the risk-free asset, \(\mu_f\) be the risk-free rate, and \(\mu_T\) be the expected return of the tangency portfolio. Then \(\omega\) solves \(0.07=\omega \mu_f+(1-\omega) \mu_T\). The following continuation of the \(\mathrm{R}\) program computes the weights for the six stocks and the risk-free asset. The last line checks that the weights sum to 1.
options(digits=3)
omega = (.07 - muP[ind]) / (mufree - muP[ind])
omega
## [1] 0.0518
1-omega
## [1] 0.948
(1-omega)*weights[ind]
## [1] -0.08645 -0.00276 0.31794 0.36382 0.30292 0.05269
omega + sum((1-omega)*weights[ind])
## [1] 1
4.3
Download monthly adjusted closing prices from Jan 2010 to Dec 2020,
for the following stocks: (Use
tseries::get.hist.quote( ..., compression=’m’) for monthly
data.)
- ENB.TO (Enbridge Inc.)
- CP.TO (Canadian Pacific Railway Ltd.)
- RCI-A.TO (Rogers Communications Inc.)
- TD.TO (The Toronto-Dominion Bank)
- L.TO (Loblaw Companies Ltd.) Calculate their net returns and assume a monthly interest rate of 0.20% throughout.
library(tseries)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
tickers = c("ENB.TO","CP.TO","RCI-A.TO","TD.TO","L.TO")
N = length(tickers)
P = vector("list", N) # list for holding prices
for (i in 1:N) {
cat("Downloading ", i, " out of ", N , "\n")
P[[i]] = get.hist.quote(instrument = tickers[i],
start = as.Date("2010-01-01"),
end=as.Date("2020-12-01"),
compression="m", quote = "AdjClose",
retclass = "zoo", quiet = T)
}
## Downloading 1 out of 5
## Downloading 2 out of 5
## Downloading 3 out of 5
## Downloading 4 out of 5
## Downloading 5 out of 5
# net returns
R = sapply(P, FUN=function(x){ as.numeric(diff(x) / timeSeries::lag(x, -1)) })
colnames(R) = tickers # assign names
head(R)
## ENB.TO CP.TO RCI-A.TO TD.TO L.TO
## [1,] 0.00496 0.00614 0.0488 0.07726 0.05158
## [2,] 0.04807 0.12699 0.0000 0.12582 0.01653
## [3,] 0.01837 0.05086 0.0592 -0.00264 0.00124
## [4,] -0.03466 -0.01670 -0.0186 -0.04195 0.04231
## [5,] 0.05019 -0.03091 -0.0137 -0.03861 -0.00797
## [6,] 0.00928 0.08071 0.0397 0.06060 0.13649
Use only the first two stocks (ENB.TO and CP.TO), and calculate the sample means and (individual) sample variances of their returns. Consider the following hypothetical values for their correlation: \(\rho=-1,-0.5,0,+0.5,+1\). For each value of \(\rho\), calculate their corresponding 2D variance-covariance matrix and plot the risk-return profiles of portfolios combining the two assets with weights \([w,(1-w)], \quad \forall w \in[-1,2]\). Plot all profile curves on the same \(\left(\mu_p, \sigma_p\right)\)-space, using a different color for each value of \(\rho\).
R2 = R[,1:2] # get first 2 stocks MU = colMeans(R2) SD = sqrt(diag(var(R2))) par(mfrow=c(1,1), mar=c(5, 4, 4, 2) + 0.1) plot(SD,MU,pch=16,cex=2, col=2, xlim=c(0,.15), ylim=c(0.005,.025)) abline(v=0, lty=2); abline(h=0, lty=2) w = seq(-2,+3,.01); W = cbind(w, 1-w) # weights MU.p = W %*% MU # portfolio means rho=c(-1, -.5, 0, .5, 1) # correlations for(i in 1:5){ COR = matrix(c(1,rho[i],rho[i],1),2,2) COV = COR * (SD%*%t(SD)) SD.p = sqrt(rowSums((W %*% COV)*W)) # portfolio st.dev. lines(SD.p, MU.p, type='l', lwd=2, col=i); } points(SD,MU,pch=16,cex=2, col=2) text(SD, MU, pch=16, colnames(R2), pos=c(1,3)) legend('topright', lwd=rep(2,5),col=1:5, c("rho=-1","rho=-.5","rho=0","rho=+.5","rho=+1") )Note that when the assets are perfectly correlated (either positively or negatively), the minimum variance portfolio has variance zero.
Consider all 5 stocks together now, and use the sample mean and sample variance-covariance matrix of their returns. Plot the efficient frontier and the capital market line on the same \(\left(\mu_p, \sigma_p\right)\)-space and report the tangency portfolio weights. (Hint: adapt the code from Example 16.6 on p. 476 of SDAFE.)
library(quadprog) COV=cov(R) MU=colMeans(R) SD=sqrt(diag(COV)) N=dim(R)[2] plot(SD, MU, pch=16, cex=1.2, col= 2, xlim=c(0,.1), ylim=c(0,.025)) abline(v=0, lty=2); abline(h=0, lty=2) text(SD, MU, tickers, cex=1, pos=4) Amat = cbind(rep(1,N),MU) mu.p = seq(-.005, .05,length=100) sd.p = mu.p; for (i in 1:length(mu.p)) { bvec=c(1,mu.p[i]) out=solve.QP(Dmat=2*COV,dvec=rep(0,N),Amat=Amat,bvec=bvec,meq=2) sd.p[i] = sqrt(out$value) } lines(sd.p,mu.p,type="l", lwd=2, col=2) # plot least variance portfolios mu.f = .002 # monthly risk-free interest rate COV.i=solve(COV) W.tang=COV.i%*%(MU-mu.f) / sum( COV.i%*%(MU-mu.f) ) # tangency ptf weights mu.tang=sum(W.tang*MU) sd.tang=sqrt(sum( (COV %*% W.tang) * W.tang ) ) points( sd.tang, mu.tang, pch=15, cex=1.3, col=2) sharpe=(mu.tang-mu.f)/sd.tang abline(mu.f,sharpe,lwd=2,lty=2,col=2)The Sharpe ratio is 0.335, and the tangency portfolio weights are:
round( t(W.tang), 4 )## ENB.TO CP.TO RCI-A.TO TD.TO L.TO ## [1,] 0.135 0.362 0.106 0.248 0.149Repeat the previous part b. (i.e. plot the efficient frontier and capital market line, and report the tangency portfolio weights) with the restriction that all weights are within the bounds: \(0 \leq w_i \leq 0.5, \forall i=1, \ldots, 5\). (Hint: adapt the code from Example 16.7 on p. 479 of SDAFE.)
plot(SD, MU, pch=16, cex=1.2, col= 2, xlim=c(0,.1), ylim=c(.0,.025)) abline(v=0, lty=2); abline(h=0, lty=2) text(SD, MU, tickers, cex=1, pos=4) Amat = cbind(rep(1,N),MU,diag(1,nrow=N),-diag(1,nrow=N)) mu.pot = seq(min(MU), max(MU),length=300) # potential mean returns mu.p = NULL # initialize portfolio standard error sd.p = NULL # initialize portfolio standard error W.p = NULL # initialize portfolio weights for (i in 1:length(mu.pot)) { bvec=c(1,mu.pot[i],rep(0,N),rep(-0.5,N)) #check whether potential mean return can be achieved with given constraints out=tryCatch( solve.QP(Dmat=2*COV,dvec=rep(0,N),Amat=Amat,bvec=bvec,meq=2), error=function(e) NULL ) #if mean return is achievable, save its st.dev. & portfolio weights if(!is.null(out)){ mu.p=c(mu.p, mu.pot[i]) sd.p=c(sd.p, sqrt(out$value)) W.p=rbind(W.p, out$solution) } } lines(sd.p, mu.p, type="l", lwd=2, col=2) # plot least variance portfolios colnames(W.p)=tickers sharpe=(mu.p-mu.f)/sd.p ind.tang=which.max(sharpe) W.tang=W.p[ind.tang,] sd.tang = sd.p[ind.tang] mu.tang = mu.p[ind.tang] points(sd.tang, mu.tang, pch=15, cex=1.3, col=2) abline(c(mu.f, sharpe[ind.tang]),lwd=2, col=2, lty=2)The set of feasible portfolios will be a subset of that of the unconstrained problem. Note that the constrained efficient frontier is not a parabola any more. However, the constrained tangency portfolio is still the same as the unconstrained one (since the unconstrained tangency portfolio weights were all within [0,.5]), and have (approximately) the same Sharpe ratio 0.335.
Note that the constrained tangency portfolio weights are (approximately) the same as the unconstrained ones:
round( W.tang, 4 )## ENB.TO CP.TO RCI-A.TO TD.TO L.TO ## 0.135 0.362 0.106 0.248 0.149
5. Factor Models
5.3 (P549, section 18.8)
Fitting Factor Models by Time Series Regression
In this section, we will start with the one-factor CAPM model of Chap. 17 and then extend this model to the three-factor Fama–French model. We will use the data set Stock_Bond_2004_to_2005.csv on the book’s website, which contains stock prices and other financial time series for the years 2004 and 2005. Data on the Fama–French factors are available at Prof. Kenneth French’s website where RF is the risk-free rate and Mkt.RF, SMB, and HML are the Fama–French factors. Go to Prof. French’s website and get the daily values of RF, Mkt.RF, SMB, and HML for the years 2004–2005. It is assumed here that you’ve put the data in a text file FamaFrenchDaily.txt. Returns on this website are expressed as percentages. Now fit the CAPM to the four stocks using the lm command. This code fits a linear regression model separately to the four responses. In each case, the independent variable is Mkt.RF.
start = as.Date("2004-01-01")
end = as.Date("2005-12-31")
stocks = read.csv("./res/datasets/Stock_Bond.csv",header=T)
stocks$Date <- as.Date(stocks$Date, format = "%d-%b-%y")
stocks_subset <- subset(stocks, Date >= start & Date <= end)
stocks_subset <- subset(stocks_subset, select = c("GM_AC", "F_AC", "UTX_AC", "MRK_AC"))
FF_data = read.table("./res/datasets/FamaFrenchDaily.txt", header = TRUE)
FF_data$date <- as.Date(as.character(FF_data$date), format = "%Y%m%d")
FF_data <- subset(FF_data, date >= start & date <= end)
FF_data = FF_data[-1, ] # delete first row since stocks_diff lost a row from differencing
stocks_diff = as.data.frame(100 * apply(log(stocks_subset),2, diff) - FF_data$RF)
names(stocks_diff) = c("GM", "Ford", "UTX", "Merck")
fit1 = lm(as.matrix(stocks_diff) ~ FF_data$Mkt.RF)
summary(fit1)
## Response GM :
##
## Call:
## lm(formula = GM ~ FF_data$Mkt.RF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.832 -0.721 0.024 0.808 15.248
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2290 0.0865 -2.65 0.0083 **
## FF_data$Mkt.RF 1.2500 0.1273 9.82 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.94 on 501 degrees of freedom
## Multiple R-squared: 0.161, Adjusted R-squared: 0.16
## F-statistic: 96.4 on 1 and 501 DF, p-value: <2e-16
##
##
## Response Ford :
##
## Call:
## lm(formula = Ford ~ FF_data$Mkt.RF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.351 -0.851 0.018 0.790 9.182
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1835 0.0676 -2.72 0.0069 **
## FF_data$Mkt.RF 1.3195 0.0995 13.26 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.51 on 501 degrees of freedom
## Multiple R-squared: 0.26, Adjusted R-squared: 0.258
## F-statistic: 176 on 1 and 501 DF, p-value: <2e-16
##
##
## Response UTX :
##
## Call:
## lm(formula = UTX ~ FF_data$Mkt.RF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.353 -0.541 -0.007 0.537 3.613
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.00219 0.03883 0.06 0.95
## FF_data$Mkt.RF 0.91932 0.05718 16.08 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.87 on 501 degrees of freedom
## Multiple R-squared: 0.34, Adjusted R-squared: 0.339
## F-statistic: 258 on 1 and 501 DF, p-value: <2e-16
##
##
## Response Merck :
##
## Call:
## lm(formula = Merck ~ FF_data$Mkt.RF)
##
## Residuals:
## Min 1Q Median 3Q Max
## -31.155 -0.481 0.059 0.709 12.307
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0888 0.0920 -0.97 0.33
## FF_data$Mkt.RF 0.6255 0.1355 4.62 5e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.06 on 501 degrees of freedom
## Multiple R-squared: 0.0408, Adjusted R-squared: 0.0389
## F-statistic: 21.3 on 1 and 501 DF, p-value: 4.99e-06
Problem 5
The CAPM predicts that all four intercepts will be zero. For each stock, using α = 0.025, can you accept the null hypothesis that its intercept is zero? Why or why not? Include the p-values with your work.
The intercept’s p-value is below 0.025 for General motors (GM) and Ford, but not for United Technologies Incorporatted (UTX) and Merck.
Of course, a p-value only shows statistical significance, not the size of an effect. However, estimated intercepts for GM and Ford are −0.23 and −0.18 and these are reasonably large in magnitude. Since they are negative, this suggests than these two stocks were overpriced; see Section 17.6.3 for discussion.
Problem 6
The CAPM also predicts that the four sets of residuals will be uncorrelated. What is the correlation matrix of the residuals? Give a 95% confidence interval for each of the six correlations. Can you accept the hypothesis that all six correlations are zero?
cor(residuals(fit1))
## GM Ford UTX Merck
## GM 1.0000 0.52016 -0.0100 -0.08776
## Ford 0.5202 1.00000 -0.0238 -0.00958
## UTX -0.0100 -0.02376 1.0000 -0.00550
## Merck -0.0878 -0.00958 -0.0055 1.00000
res = residuals(fit1)
cor.test(res[,"GM"],res[,"Ford"])
##
## Pearson's product-moment correlation
##
## data: res[, "GM"] and res[, "Ford"]
## t = 14, df = 501, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.453 0.581
## sample estimates:
## cor
## 0.52
cor.test(res[,"GM"],res[,"Merck"])
##
## Pearson's product-moment correlation
##
## data: res[, "GM"] and res[, "Merck"]
## t = -2, df = 501, p-value = 0.05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.17385 -0.00033
## sample estimates:
## cor
## -0.0878
All correlations are reasonably close to 0 (less than 0.1 in magnitude) except the correlation between GM and Ford residuals. That correlation is 0.52 and has a very small p-value.
The correlation between GM and Merck residuals is −0.0878 and is statistically significant at 0.05 but might be too small to be of practical significance.
Problem 7
Regardless of your answer to Problem 6, assume for now that the residuals are uncorrelated. Then use the CAPM to estimate the covariance matrix of the excess returns on the four stocks. Compare this estimate with the sample covariance matrix of the excess returns. Do you see any large discrepancies between the two estimates of the covariance matrix?
attach(FF_data)
sigF = var(Mkt.RF)
bbeta = as.matrix(fit1$coef)
bbeta = bbeta[-1,] # delete intercepts so bbeta has the four slopes
n=dim(stocks_diff)[1]
sigeps = as.matrix((var(as.matrix(res))))
sigeps = diag(as.matrix(sigeps))
sigeps = diag(sigeps,nrow=4)
cov_equities = sigF* bbeta %*% t(bbeta) + sigeps
cov_equities
## GM Ford UTX Merck
## [1,] 4.464 0.761 0.530 0.360
## [2,] 0.761 3.090 0.559 0.381
## [3,] 0.530 0.559 1.145 0.265
## [4,] 0.360 0.381 0.265 4.423
cov(stocks_diff)
## GM Ford UTX Merck
## GM 4.4641 2.282 0.513 0.0108
## Ford 2.2825 3.090 0.528 0.3507
## UTX 0.5130 0.528 1.145 0.2553
## Merck 0.0108 0.351 0.255 4.4228
We see that the estimated covariance matrix using the CAPM is similar to the sample covariance matrix, with the exception of the covariance between GM and Ford. Since these two stocks have a high residual correlation and the CAPM assumes that the residual correlation is 0, it is not surprising that the CAPM estimated covariance matrix severely underestimates the correlation between GM and Ford.
Problem 8
Next, you will fit the Fama–French three-factor model. Run the following R code, which is much like the previous code except that the regression model has two additional predictor variables, SMB and HML.
fit2 = lm(as.matrix(stocks_diff) ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
summary(fit2)
## Response GM :
##
## Call:
## lm(formula = GM ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.759 -0.727 -0.033 0.775 15.007
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2511 0.0865 -2.90 0.0038 **
## FF_data$Mkt.RF 1.3889 0.1517 9.16 <2e-16 ***
## FF_data$SMB -0.2504 0.2236 -1.12 0.2633
## FF_data$HML 0.6006 0.2643 2.27 0.0235 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.93 on 499 degrees of freedom
## Multiple R-squared: 0.173, Adjusted R-squared: 0.168
## F-statistic: 34.9 on 3 and 499 DF, p-value: <2e-16
##
##
## Response Ford :
##
## Call:
## lm(formula = Ford ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.344 -0.832 -0.010 0.793 9.198
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.1951 0.0679 -2.87 0.0042 **
## FF_data$Mkt.RF 1.3511 0.1191 11.35 <2e-16 ***
## FF_data$SMB -0.0157 0.1756 -0.09 0.9288
## FF_data$HML 0.3412 0.2074 1.64 0.1006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.51 on 499 degrees of freedom
## Multiple R-squared: 0.264, Adjusted R-squared: 0.26
## F-statistic: 59.7 on 3 and 499 DF, p-value: <2e-16
##
##
## Response UTX :
##
## Call:
## lm(formula = UTX ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.269 -0.537 0.008 0.493 3.429
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.99e-06 3.88e-02 0.00 0.9999
## FF_data$Mkt.RF 1.03e+00 6.80e-02 15.12 <2e-16 ***
## FF_data$SMB -2.93e-01 1.00e-01 -2.92 0.0037 **
## FF_data$HML -9.59e-04 1.19e-01 -0.01 0.9935
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.864 on 499 degrees of freedom
## Multiple R-squared: 0.352, Adjusted R-squared: 0.348
## F-statistic: 90.2 on 3 and 499 DF, p-value: <2e-16
##
##
## Response Merck :
##
## Call:
## lm(formula = Merck ~ FF_data$Mkt.RF + FF_data$SMB + FF_data$HML)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.627 -0.454 0.106 0.691 12.310
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0598 0.0915 -0.65 0.51337
## FF_data$Mkt.RF 0.7093 0.1605 4.42 1.2e-05 ***
## FF_data$SMB -0.4174 0.2366 -1.76 0.07832 .
## FF_data$HML -0.9559 0.2796 -3.42 0.00068 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.04 on 499 degrees of freedom
## Multiple R-squared: 0.066, Adjusted R-squared: 0.0604
## F-statistic: 11.8 on 3 and 499 DF, p-value: 1.86e-07
The CAPM predicts that for each stock, the slope (beta) for SMB and HML will be zero. Explain why the CAPM makes this prediction. Do you accept this null hypothesis? Why or why not?
Problem 9
If the Fama–French model explains all covariances between the returns, then the correlation matrix of the residuals should be diagonal. What is the estimated correlations matrix? Would you accept the hypothesis that the correlations are all zero?
Problem 10
Which model, CAPM or Fama–French, has the smaller value of AIC? Which has the smaller value of BIC? What do you conclude from this?
Problem 11
What is the covariance matrix of the three Fama–French factors?
Problem 12
In this problem, Stocks 1 and 2 are two stocks, not necessarily in the Stock_FX_Bond_2004_to_2005.csv data set. Suppose that Stock 1 has betas of 0.5, 0.4, and −0.1 with respect to the three factors in the Fama–French model and a residual variance of 23.0. Suppose also that Stock 2 has betas of 0.6, 0.15, and 0.7 with respect to the three factors and a residual variance of 37.0. Regardless of your answer to Problem 9, when doing this problem, assume that the three factors do account for all covariances.
Use the Fama–French model to estimate the variance of the excess return on Stock 1.
Use the Fama–French model to estimate the variance of the excess return on Stock 2.
Use the Fama–French model to estimate the covariance between the excess returns on Stock 1 and Stock 2.
Statistical Factor Models
This section applies statistical factor analysis to the log returns of 10 stocks in the data set Stock_FX_Bond.csv. The data set contains adjusted closing (AC) prices of the stocks, as well as daily volumes and other information that we will not use here. The following R code will read the data, compute the log returns, and fit a two-factor model. Note that factanal works with the correlation matrix or, equivalently, with standardized variables.
dat = read.csv("./res/datasets/Stock_FX_Bond.csv")
stocks_ac = dat[ , c(3, 5, 7, 9, 11, 13, 15, 17)]
n = length(stocks_ac[ , 1])
stocks_returns = log(stocks_ac[-1, ] / stocks_ac[-n, ])
fact = factanal(stocks_returns, factors = 2, rotation = "none")
print(fact)
##
## Call:
## factanal(x = stocks_returns, factors = 2, rotation = "none")
##
## Uniquenesses:
## GM_AC F_AC UTX_AC CAT_AC MRK_AC PFE_AC IBM_AC MSFT_AC
## 0.399 0.423 0.718 0.714 0.519 0.410 0.760 0.749
##
## Loadings:
## Factor1 Factor2
## GM_AC 0.693 -0.348
## F_AC 0.692 -0.313
## UTX_AC 0.531
## CAT_AC 0.529
## MRK_AC 0.551 0.421
## PFE_AC 0.574 0.511
## IBM_AC 0.490
## MSFT_AC 0.499
##
## Factor1 Factor2
## SS loadings 2.64 0.666
## Proportion Var 0.33 0.083
## Cumulative Var 0.33 0.414
##
## Test of the hypothesis that 2 factors are sufficient.
## The chi square statistic is 565 on 13 degrees of freedom.
## The p-value is 2.6e-112
Loadings less than the parameter cutoff are not printed. The default value of cutoff is 0.1, but you can change it as in “print(fact,cutoff = 0.01)” or “print(fact, cutoff = 0)”.
Problem 13
What are the factor loadings? What are the variances of the unique risks for Ford and General Motors?
We see from the output that the uniquenesses for Ford and GM are 0.423 and 0.399, respectively. The loadings are printed below (also shown in the output).
unique = as.numeric(fact$unique)
loadings = matrix(as.numeric(loadings(fact)), ncol = 2)
loadings
## [,1] [,2]
## [1,] 0.693 -0.34769
## [2,] 0.692 -0.31339
## [3,] 0.531 -0.02819
## [4,] 0.529 -0.07850
## [5,] 0.551 0.42090
## [6,] 0.574 0.51112
## [7,] 0.490 0.00424
## [8,] 0.499 0.04400
Problem 14
Does the likelihood ratio test suggest that two factors are enough? If not, what is the minimum number of factors that seems sufficient? The following code will extract the loadings and uniquenesses.
# try increasing the number of factors
factanal(stocks_returns, factors = 3, rotation = "none")
##
## Call:
## factanal(x = stocks_returns, factors = 3, rotation = "none")
##
## Uniquenesses:
## GM_AC F_AC UTX_AC CAT_AC MRK_AC PFE_AC IBM_AC MSFT_AC
## 0.402 0.368 0.713 0.703 0.510 0.386 0.577 0.622
##
## Loadings:
## Factor1 Factor2 Factor3
## GM_AC 0.688 -0.339
## F_AC 0.704 -0.334 -0.160
## UTX_AC 0.528
## CAT_AC 0.528 0.112
## MRK_AC 0.546 0.420 -0.126
## PFE_AC 0.572 0.524 -0.111
## IBM_AC 0.518 0.392
## MSFT_AC 0.520 0.324
##
## Factor1 Factor2 Factor3
## SS loadings 2.689 0.687 0.343
## Proportion Var 0.336 0.086 0.043
## Cumulative Var 0.336 0.422 0.465
##
## Test of the hypothesis that 3 factors are sufficient.
## The chi square statistic is 162 on 7 degrees of freedom.
## The p-value is 1.06e-31
factanal(stocks_returns, factors = 4, rotation = "none")
##
## Call:
## factanal(x = stocks_returns, factors = 4, rotation = "none")
##
## Uniquenesses:
## GM_AC F_AC UTX_AC CAT_AC MRK_AC PFE_AC IBM_AC MSFT_AC
## 0.411 0.349 0.077 0.702 0.496 0.402 0.540 0.620
##
## Loadings:
## Factor1 Factor2 Factor3 Factor4
## GM_AC 0.489 0.484 -0.333
## F_AC 0.496 0.516 -0.345 -0.141
## UTX_AC 0.936 -0.217
## CAT_AC 0.486 0.229
## MRK_AC 0.373 0.406 0.428 -0.129
## PFE_AC 0.392 0.413 0.513 -0.107
## IBM_AC 0.390 0.343 0.435
## MSFT_AC 0.379 0.351 0.329
##
## Factor1 Factor2 Factor3 Factor4
## SS loadings 2.185 1.176 0.685 0.356
## Proportion Var 0.273 0.147 0.086 0.044
## Cumulative Var 0.273 0.420 0.506 0.550
##
## Test of the hypothesis that 4 factors are sufficient.
## The chi square statistic is 0.3 on 2 degrees of freedom.
## The p-value is 0.86
The results of the likelihood ratio tests below strongly suggest that there are more than two factors. It seems that four factors are sufficient, but not 3.
Problem 15
Regardless of your answer to Problem 14, use the two-factor model to estimate the correlation of the log returns for Ford and IBM
loadings = matrix(as.numeric(loadings(fact)),ncol=2)
unique = as.numeric(fact$unique)
options(digits=2)
loadings %*% t(loadings) + diag(unique)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
## [1,] 1.00 0.59 0.38 0.39 0.24 0.22 0.34 0.33
## [2,] 0.59 1.00 0.38 0.39 0.25 0.24 0.34 0.33
## [3,] 0.38 0.38 1.00 0.28 0.28 0.29 0.26 0.26
## [4,] 0.39 0.39 0.28 1.00 0.26 0.26 0.26 0.26
## [5,] 0.24 0.25 0.28 0.26 1.00 0.53 0.27 0.29
## [6,] 0.22 0.24 0.29 0.26 0.53 1.00 0.28 0.31
## [7,] 0.34 0.34 0.26 0.26 0.27 0.28 1.00 0.24
## [8,] 0.33 0.33 0.26 0.26 0.29 0.31 0.24 1.00
cor(stocks_returns)
## GM_AC F_AC UTX_AC CAT_AC MRK_AC PFE_AC IBM_AC MSFT_AC
## GM_AC 1.00 0.62 0.35 0.36 0.25 0.23 0.32 0.31
## F_AC 0.62 1.00 0.35 0.37 0.26 0.25 0.30 0.30
## UTX_AC 0.35 0.35 1.00 0.40 0.26 0.28 0.29 0.28
## CAT_AC 0.36 0.37 0.40 1.00 0.24 0.25 0.30 0.29
## MRK_AC 0.25 0.26 0.26 0.24 1.00 0.55 0.24 0.27
## PFE_AC 0.23 0.25 0.28 0.25 0.55 1.00 0.26 0.29
## IBM_AC 0.32 0.30 0.29 0.30 0.24 0.26 1.00 0.41
## MSFT_AC 0.31 0.30 0.28 0.29 0.27 0.29 0.41 1.00
5.4
For this problem you will use regression to identify the composition of various mutual funds.
(a)
Download the adjusted daily closing prices from Jan 1 2020 to Dec 31
2022 for the 5 mutual funds below (use
tseries::get.hist.quote() for each ticker):
• FCNTX: Fidelity Contrafund
• PIODX: Pioneer A
• AIVSX: American Funds Invmt Co of Amer A
• PRBLX: Parnassus Core Equity Investor
• VFIAX: Vanguard 500 Index Admiral
Note that each of these funds has at least 90% of their weight in the US stocks market. You can actually check the composition of the investment over different stock sectors from Yahoo Finance, under the fund’s holdings tab; e.g. for FCNTX at https://finance.yahoo.com/quote/FCNTX/holdings.
library(zoo)
library(tseries)
MF.names=c('FCNTX','PIODX','AIVSX','PRBLX','VFIAX')
N.MF=length(MF.names)
S=list()
for(i in 1:N.MF){
S[[i]]=get.hist.quote(MF.names[i], start='2022-01-01',
end='2022-12-31', quote='AdjClose', quiet = TRUE)
}
R=lapply(S, FUN = function(x){ diff(x) / timeSeries::lag(x,-1) }) # calculate MF returns
RY=matrix(unlist(R),ncol=N.MF) # bind returns in a matrix
colnames(RY)=MF.names
(b)
Assume you do not have any information about the investment strategy
of the funds.Download the daily prices and calculate returns of the
following EFTs, which track different sectors of the economy:
• XLB: Basic Materials
• XLY: Consumer Cyclical
• XLF: Financial Services
• VNQ: Real Estate
• XLP: Consumer Defensive
• XLV: Healthcare
• XLU: Utilities
• XTL: Communication Services
• XLE: Energy
• XLI: Industrials
• XLK: Technology
Regress each of the mutual fund returns on the above ETF returns and create barplots of the estimated beta coefficients. Do these accurately reflect the allocation over the different sectors (as described in Yahoo Finance)?
ETF.names=c('XLB','XLY','XLF','VNQ','XLP','XLV','XLU','XTL','XLE','XLI','XLK')
sectors = c("Basic Materials","Consumer Cyclical","Financial Services",
"Real Estate","Consumer Defensive", "Healthcare", "Utilities",
"Communication Services", "Energy", "Industrials", "Technology")
N.ETF=length(ETF.names)
S=list()
for(i in 1:N.ETF){
S[[i]]=get.hist.quote(ETF.names[i], start='2022-01-01',
end='2022-12-31', quote='AdjClose', quiet = TRUE)
}
## Warning: XLE download failed; trying again.
R=lapply(S, FUN = function(x){ diff(x) / timeSeries::lag(x,-1) }) # calculate ETF returns
RX=matrix(unlist(R),ncol=N.ETF) # bind returns in a matrix
colnames(RX)=ETF.names
out=list();
for(i in 1:N.MF){
out[[i]]=lm(RY[,i] ~ RX) # regress mutual funds returns on ETF returns
weights = out[[i]]$coef[-1]
barplot( rev(weights), names.arg = rev(sectors), main = MF.names[i],
horiz = TRUE, las = 2, cex.names = .8)
}
The barplots of the regression coefficients (betas) roughly follow the sector weightings for each fund. Nevertheless, they are not always close in actual value (e.g. in some cases the betas are negative, even though weightings are positive). The differences can be due to the fact that we use ETFs as proxies for a sector, but the actual holding of the fund within the sector might be different. Moreover, there will be estimation error in our regression model, which is only based on the last year’s returns.
Nevertheless, it is quite impressive that we can (approximately) identify the strategy of a fund, without knowing anything beyond its past returns. This approach works because of the linear formula for net portfolio returns: \[ R_p = w_1 R_1 + \cdots + w_N R_N\] Regressing portfolio returns on other assets, we can estimate the weights, assuming the portfolio composition is constant.
(c)
Compare the performance of the mutual funds to that of a portfolio of ETFs by reporting the value of Jensen’s alpha (based on the regressions from the previous part) and its corresponding p-value.
alpha=p.val=rep(0,N.MF)
for(i in 1:N.MF){
alpha[i]=out[[i]]$coefficients[1]*250 # annualized Jensen alpha
p.val[i]=summary(out[[i]])$coefficients[1,4]
}
cbind(alpha, p.val)
## alpha p.val
## [1,] -0.103 0.078
## [2,] -0.034 0.370
## [3,] -0.039 0.250
## [4,] -0.025 0.468
## [5,] -0.016 0.348
All the funds’ alphas are negative, although their p-values are not very small. A likely cause for this is that funds charge a fee which consistently eats up some of the returns of their constituent assets. ETFs typically have lower fees than mutual funds, but our regression does not account for transaction costs (it is more costly to buy multiple assets than a single one), so the comparison is more nuanced.
Note that you can find Jensen alphas and other performance measures (e.g., Sharpe & Treynor Ratios) for assets in Yahoo! Finance under the risk tab. These metrics are based on the CAPM/Market factor model, by regressing the asset’s returns on a proxy for the market return (e.g., S&P500).
6. Risk Management
6.7 Univariate VaR and ES
In this section we will compare VaR and ES parametric (unconditional) estimates with those from using ARMA+GARCH (conditional) models. Consider the daily returns for Coca-Cola stock from January 2007 to November 2012.
CokePepsi = read.table("./res/datasets/CokePepsi.csv", header=T)
price = CokePepsi[,1]
returns = diff(price)/lag(price)[-1]
ts.plot(returns)
First, assume that the returns are iid and follow a t-distribution. Run the following commands to get parameter estimates in R.
S = 4000
alpha = 0.05
res = fitdistr(returns,'t')
## Warning in dt((x - m)/s, df, log = TRUE): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in dt((x - m)/s, df, log = TRUE): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
## Warning in log(s): NaNs produced
mu = res$estimate['m']
lambda = res$estimate['s']
nu = res$estimate['df']
qt(alpha, df=nu)
## [1] -2.3
dt(qt(alpha, df=nu), df=nu)
## [1] 0.048
For an investment of $4,000, what are estimates of \(VaR^t(0.05)\) and \(ES^t(0.05)\)?
mu = as.numeric(res$estimate['m'])
lambda = as.numeric(res$estimate['s'])
nu = as.numeric(res$estimate['df'])
qt(alpha, df=nu)
## [1] -2.3
dt(qt(alpha, df=nu), df=nu)
## [1] 0.048
Finv = mu + lambda * qt(alpha, df=nu)
VaR = -S * Finv
options(digits=4)
VaR
## [1] 75.31
den = dt(qt(alpha, df=nu), df=nu)
ES = S * (-mu + lambda*(den/alpha) * (nu+qt(alpha, df=nu)^2 )/(nu-1))
ES
## [1] 122.1
(Extra) Fit a ARMA(0,0)+GARCH(1,1) model to the returns and calculate one step forecasts.
# library(fGarch) # for qstd() function
# library(rugarch)
# garch.t = ugarchspec(mean.model=list(arma0rder=c(0,0)) ,
# variance.model=list(garchOrder=c (1,1)),
# distribution.model = "std")
# KO.garch.t = ugarchfit(data=returns, spec=garch.t)
# show(KO.garch.t)
# plot(KO.garch.t, which = 2)
# pred = ugarchforecast(KO.garch.t, data-returns, n.ahead=1); pred
# fitted(pred); sigma(pred)
# nu = as.numeric(coef(KO.garch.t) [5])
# q = qstd(alpha, mean = fitted(pred), sd = sigma(pred), nu = nu) ; q
# sigma(pred) /sqrt( (nu) / (nu-2) )
# qt(alpha, df=nu)
# dt(qt(alpha, df=nu) , df=nu)
Midterms
2021 Q2 (Generalized Pareto Distribution)
For this question you will generate GPD random variates and verify that their conditional excess distribution is also GPD.
(a) [10 points]
Simulate \(n=10,000\) random variates from a \(\operatorname{GPD}(\gamma=.5, \sigma=1)\) distribution using the inverse CDF method from Q1.(d), and plot the histogram of their log-values.?
set.seed(123)
n = 10000
gamma = .5
sigma = 1
X = sigma/gamma * ( (1-runif(n))^(-gamma) - 1 )
hist( log( X ) )
(b) [15 points]
Simulate another 10, 000 random variates from \(\operatorname{GPD}(\gamma=.5, \sigma=1)\), but now use the mixture method from Q1.(e). Create a QQ-plot of the two sets of variates (from this and the previous part0), and comment whether they seem to come from the same distribution? (Hint: Use qqplot \((\mathrm{X}, \mathrm{Y})\), where \(\mathrm{X}, \mathrm{Y}\) are the values you generated.)
alpha = 1 / gamma
beta = sigma / gamma
L = rgamma( n, shape = alpha, rate = beta )
Y = rexp( n, rate = L)
qqplot(X,Y); abline( c(0,1) )
(c) [15 points]
Take the values from Q2.(a) and calculate their exceedances above 10, i.e. \(Z=X-10\) for values of \(X\) that are greater than 10. From Q1.(b), we know that \(Z\) must follow \(G P D(\gamma, \sigma+10 \gamma)\). Verify that the exceedances follow this distribution, by creating a QQ-plot of \(Z\)-values versus their theoretical quantiles.
u = 10
# theoretical quantiles
Z = (X-u)[X>u] # values of X-u that satisfy X>u
Zs = sort(Z)
# sample quantiles
nZ = length(Z)
sigma_u = sigma + gamma * u
Q = sigma_u/gamma * ( (1 - ppoints(nZ) )^(-gamma) - 1 )
# ppoints(n) creates n increments of size 1/n
plot( Zs, Q ); abline(0,1)
(d) [10 points]
Finally, verify through simulation that conditional exceedances of fat tail distributions with tail index \(\alpha\) approach a GPD with \(\gamma=1 / \alpha\). Generate \(n=10,000\) values from a \(t(d f=4)\) distribution (i.e. \(\alpha=4\) ) and take their absolute value (to avoid wasting values by symmetry). Calculate the exceedances above 10 again, and create their QQ-plot versus quantiles from a \(\operatorname{GPD}(\gamma=1 / \alpha, \sigma=1)\). (Note: the points in your QQ-plot should lie close to a straight line, but the slope does not need to be equal to 1 because of the arbitrary GPD scale \(\sigma=1\) ).
# theoretical quantiles
W = abs( rt( 100000, 4 ) )
Ws = (W-u)[W>u]
Ws = sort(Ws)
# sample quantiles
nW = length(W)
sigma = 1
gamma = 1/4
Q = (1+u) * ( (1 - ppoints(nW))^(-gamma) - 1 ) # why (1+u)?
# Q = ((sigma + gamma * u) / gamma) * ( (1 - ppoints(nW))^(-gamma) - 1 )
qqplot( Ws, Q ); abline(0, 1)
2020 Q3 (Cauchy distribution)
You have to perform a simulation experiment to confirm that the t(1) or Cauchy distribution is stable, but not all heavy-tailed distributions are.
part (a)
[4 points] Simulate 10, 000 Cauchy random variates and plot their cumulative mean versus the sample size, as in the plot below
n = 10000
set.seed(123)
x = rcauchy(n)
m = cumsum(x)/(1:n)
plot(m, type = 'l' )
part (b)
[8 points] Repeat the following experiment N = 1000 times: simulate n = 100 Cauchy variates and calculate their mean (in the end you should have 1,000).
From question 1, we know these means should be also Cauchy distributed. Verify this by creating a QQ-plot of the simulated values versus their theoretical quantiles, and comment on the fit.
Hint: Use qcauchy and ppoints() to generate theoretical quantiles from the Cauchy distribution.)
N = 1000; n = 100
x = replicate( n = N, expr = mean(rcauchy(n = n)) )
qqplot( sort(x), qcauchy( p = ppoints(N)) ); abline(c(0,1))
part (c)
[8 points] Now consider the Pareto distribution, with CDF \(F_X(x)=1-(1 / x)^\alpha, \forall x>1\) Use the inverse CDF method to generate 10,000 Pareto random variates with tail index \(\alpha=3\), and plot their cumulative mean vs sample size, as in part (a).
n = 10000
alpha = 3
u = runif(n)
y = (1 - u)^(-1/alpha)
my = cumsum(y)/(1:n)
plot(my, type = 'l' )
abline(h=1.5, col = 2, lty=2)
2020 Q4 (Factor Model)
Consider the 3-factor model for the returns of two assets: \[ \left[\begin{array}{l} R_1 \\ R_2 \end{array}\right]=\left[\begin{array}{l} .03 \\ .04 \end{array}\right]+\left[\begin{array}{lll} .1 & .2 & .3 \\ .3 & .2 & .1 \end{array}\right]\left[\begin{array}{l} F_1 \\ F_2 \\ F_3 \end{array}\right]+\left[\begin{array}{l} e_1 \\ e_2 \end{array}\right] \Leftrightarrow \boldsymbol{R}=\boldsymbol{\mu}+\boldsymbol{\beta}^{\top} \boldsymbol{F}+\boldsymbol{e} \] where \(\boldsymbol{F} \sim N_3(\mathbf{0}, \boldsymbol{I})\) and \(\boldsymbol{e} \sim N_2\left(\mathbf{0}, \boldsymbol{\Sigma}_e\right)\) with \(\boldsymbol{\Sigma}_e=\left[\begin{array}{cc}.1 & 0 \\ 0 & .2\end{array}\right]\) and \(\boldsymbol{F} \perp \boldsymbol{e}\).
part (a)
[5 points] Calculate the variance-covariance matrix of the returns \(\boldsymbol{\Sigma}_R=\mathbb{V}[\boldsymbol{R}]\).
beta = matrix( c(1,2,3,3,2,1)/10, nrow = 2, byrow = TRUE )
S_e = diag( c(.1,.2) )
S_R = beta %*% t(beta) + S_e
part (b)
[5 points] Calculate the mean and variance of the return of the minimum-variance portfolio.
Si = solve(S_R) # Sigma inverse
w = rowSums(Si)/sum(Si)
mu = c(.03, .04)
mvp_mu = sum(w*mu)
mvp_s = t(w) %*% S_R %*% w # also equal to 1/sum(Si)
part (c)
[10 points] Assume the risk-free rate is .015, and find the mean and variance of the tangency portfolio.
rf = .015
w_tp = Si %*% (mu - rf) / sum( Si %*% (mu - rf) )
sum( w_tp * mu ) # tangency portfolio mean
## [1] 0.03634
t(w_tp) %*% S_R %*% w_tp # tangency portfolio variance
## [,1]
## [1,] 0.2152
part (d)
[10 points] Assume you could eliminate one of the factors (i.e. set \(\mathbb{V}\left[F_i\right]=0\) ). Which factor would you choose to eliminate in order to further reduce the variance of the minimumvariance portfolio? (justify your answer with a calculation & comparison)
mvp_vars = rep(0,3)
for(i in 1:3){
beta_t = beta[,-i] # t for two factor
S_t = beta_t %*% t(beta_t) + S_e
Si_t = solve(S_t)
mvp_vars[i] = 1/sum(Si_t) # min variance
}
which.min(mvp_vars)
## [1] 3
7. Betting Strategies
7.3 (Kelly for predictable returns)
The “Efficient Markets Hypothesis” presumes that returns are unpredictable. For this question, however, you will assume that returns are predictable and they actually follow an AR(1) process \(\left(R_t-\mu\right)=\varphi\left(R_{t-1}-\mu\right)+\varepsilon_t\), where \(\varepsilon_t \sim^{i i d} N\left(0, \sigma_{\varepsilon}^2\right)\). The model dictates that the conditional mean of the next return is a linear function of the current return, i.e. \(\mathbb{E}\left[R_t-\mu \mid R_{t-1}\right]=\varphi\left(R_{t-1}-\mu\right)\). Based on this observation, one can modify the Kelly criterion for investing that we saw in class to allow for dynamically changing fractions, based on the conditional mean & variance of returns at each point. Remember that the Kelly-optimal fraction for i.i.d. returns with mean \(\mu\), variance \(\sigma^2\) and constant risk-free rate \(r_f\) is \(f^* \approx\left(1+r_f\right) \frac{\mu-r_f}{\sigma^2}\). In the dynamic setting, substitute \(\mu\) and \(\sigma^2\) with the conditional mean and variance of the return, i.e. \[ f_{t-1}^* \approx\left(1+r_f\right) \frac{\mathbb{E}\left[R_t \mid R_{t-1}\right]-r_f}{\mathbb{V}\left[R_t \mid R_{t-1}\right]} \]
(a)
Simulate \(n=1,000\) paths of \(m=252\) daily returns each, from an AR(1) model with \(R_0=0, \varphi=.25, \mu=.05 / 252, \sigma_{\varepsilon}=.2 / \sqrt{252}, r_f=.02 / 252\).
n=1000; m=252
phi=.25; mu=.05/252; sigma=.4/sqrt(252); rf=.02/252
V0=1
set.seed(1234567890)
R=matrix(0,m,n)
Z=matrix(rnorm(n*m),m,n)
# Initialize returns at stationary distribution
R[1,]=mu+(0-mu)*phi+sigma*Z[1,]
for(i in 2:m){
R[i,]=mu+(R[i-1,]-mu)*phi+sigma*Z[i,]
}
# To avoid the loop (& speed up computation) you can use the filter() function:
# R = apply( mu*(1-phi) + sigma*Z, 2, filter, filter=phi, method="recursive")
(b)
Apply the static Kelly criterion based on the stationary mean \(\mathbb{E}\left[R_t\right]=\mu\) and variance \(\mathbb{V}\left[R_t\right]=\frac{\sigma_{\varepsilon}^2}{\left(1-\varphi^2\right)}\) of the \(\mathrm{AR}(1)\) process, and calculate the resulting wealth process for each path. Plot the wealth of the strategy on the first path, and the histogram of the final log-wealths across all paths.
sigma.b=sigma/sqrt(1-phi^2)
(f.b=(1+rf)*(mu-rf)/(sigma.b^2)) # static Kelly fraction
## [1] 0.1758
V.b=matrix(0,m,n)
V.b[1,]=V0*(1+rf+f.b*(R[1,]-rf))
for(i in 2:m){
V.b[i,]=V.b[i-1,]*(1+rf+f.b*(R[i,]-rf))
}
# To avoid the loop use
# V.b=apply( 1+rf+f.b*(R-rf), 2, cumprod)
plot(V.b[,1], type='l')
hist( V.b[m,], main = "Static Kelly", xlab="Vn" )
Using the “static” Kelly criterion, we do not really take advantage of the predictability of the process, just of the fact that the asset returns are on average better than the risk-free rate. The average final wealth over 1 year is a modest 1.0247.
(c)
Apply the dynamic Kelly criterion and calculate the resulting wealth process for each path. Plot the first path and the histogram of the final log-wealths across all paths, and compare them to those of the previous part.
V.c=matrix(1,m,n)
# Dynamic Kelly fraction
mu_i=(0-mu)*phi
(f.c_i=(1+rf)*mu_i/(sigma^2)) # if f<0 -> short asset
## [1] -0.07813
V.c[1,]=V0*(1+rf+f.c_i*(R[1,]-rf))
for(i in 2:m){
# Dynamic Kelly fraction
mu_i=(R[i-1,]-mu)*phi
f.c_i=(1+rf)*mu_i/(sigma^2) # changing fraction
V.c[i,]=V.c[i-1,]*(1+rf+f.c_i*(R[i,]-rf))
}
# To avoid the loop use
# mu_i=( rbind(rep(0,n),R[-m,])-mu)*phi
# f.c_i=(1+rf)*mu_i/(sigma^2+mu_i^2)
# V.c=apply( 1+rf+f.c_i*(R-rf), 2, cumprod)
# Plot wealth of 1st path
plot(V.c[,1], type='l', main = "Wealth path")
lines(V.b[,1], type='l', col=2)
legend( "topleft", lwd = 2, col = 1:2, c("dynamic","static"))
# Plot log-wealth, for better comparison.
plot(V.c[,1], type='l', log = "y", main = "log-Wealth path")
lines(V.b[,1], type='l', col=2)
legend( "topleft", lwd = 2, col = 1:2, c("dynamic","static"))
The dynamic Kelly wealth is growing much faster than the static one. That’s because it can take advantage of the predictability of the price on each day, and even when the price goes down, by using a “negative” fraction (i.e. shorting the asset). However, it is possible for the dynamic Kelly wealth to become negative if the optimal fraction is >1 (i.e., you use leverage). In fact, 22.4% of the dynamic Kelly paths had negative wealth (i.e. ruin/bankruptcy) at some point. But for those paths that did not go negative, the dynamic strategy led to much higher terminal wealth.
par(mfrow=c(1,2))
hist(log(V.b[m,]), main = "Static Kelly", xlab = "log(Vn)" )
hist(log(V.c[m, V.c[m,]>0]), main = "Dynamic Kelly", xlab = "log(Vn[Vn>0])" )
One can moderate and improve the dynamic strategy by constraining the invested fraction to be always \(|f^*|<1\) (i.e. avoid leverage). This technically removes the risk of ruin while still giving superior results:
V.c.constr=matrix(1,m,n)
mu_i=(0-mu)*phi
f.c_i= min( max( (1+rf)*mu_i/(sigma^2), -1 ), 1)
V.c[1,]=V0*(1+rf+f.c_i*(R[1,]-rf))
for(i in 2:m){
# Dynamic Kelly fraction
mu_i=(R[i-1,]-mu)*phi
f.c_i= pmax( pmin( (1+rf)*mu_i/(sigma^2), 1) , -1)
V.c.constr[i,]=V.c.constr[i-1,]*(1+rf+f.c_i*(R[i,]-rf))
}
par(mfrow=c(1,2))
hist(log(V.b[m,]), main = "Static Kelly", xlab = "log(Vn)" )
hist(log(V.c.constr[m,]), main = "Dynamic Kelly", xlab = "log(Vn)" )
8. Statistical Arbitrage
When is \(x/y - 1\) is a good approximation for \(\log x - \log y\)?
f1 = function(x, y){x/y - 1}
f2 = function(x, y){log(x) - log(y)}
n = 100
x <- y <- seq(0.5, 100.5, length = n)
z1 = outer(x, y, f1) # 100 x 100 matrix
z2 = outer(x, y, f2)
z3 = outer(x, y, function(x, y) x)
# Plot the first surface
persp(x, y, z1, theta = 30, phi = 30, col = "lightblue", shade = 0.5, border = NA,
xlab = "X", ylab = "Y", zlab = "Z")
# Add the second surface to the same plot
par(new = TRUE)
persp(x, y, z2, theta = 30, phi = 30, col = "lightgreen", shade = 0.5, border = NA,
xlab = "", ylab = "", zlab = "", axes = FALSE, box = FALSE)
# Add the third surface (y=x) to the same plot
par(new = TRUE)
persp(x, y, z3, theta = 30, phi = 30, col = "red", shade = 0.5, border = NA,
xlab = "", ylab = "", zlab = "", axes = FALSE, box = FALSE)
legend("topright", c("x/y - 1", "log(x/y)", "y=x"), fill = c("lightblue", "lightgreen", "red"))
# write function with param specifying offset (for a square matrix)
diag_offset <- function(mat, offset){
n = dim(mat)[1]
if (offset >= 0){
m = mat[1:(n-offset), (1+offset):n]
}
else {
m = mat[(1+offset):n, 1:(n-offset)]
}
diag(m)
}
offset = 1 # range: -(n-1) to (n-1)
# Create plot
plot(diag_offset(z1, offset) - diag_offset(z2, offset), type = 'l')
for (i in -(100-1):(100-1)) {
lines(diag_offset(z1, offset) - diag_offset(z2, offset), col = abs(i), lwd = 2)
}
# Load required library
library(ggplot2)
library(reshape2)
# Create the heatmap
m = z1-z2
ggplot(data = melt(m), aes(x = Var2, y = Var1, fill = value)) +
geom_tile() +
scale_fill_gradient(low = 'blue', high = "red") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
# Shows y cannot be too small
# Take subset to exclude small y values, replot
m = m[1:100, 10:100]
ggplot(data = melt(m), aes(x = Var2, y = Var1, fill = value)) +
geom_tile() +
scale_fill_gradient(low = 'blue', high = "red") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
# we see that the difference between z1 and z2 is lower when x is close to y
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:MASS':
##
## select
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
# Define x, y ranges
x <- seq(0.01, 10, length.out = 100)
y <- seq(0.01, 10, length.out = 100)
# Define functions
f1 <- function(x, y) x/(y+1)
f2 <- function(x, y) log(x) - log(y)
f3 <- function(x, y) x - y
# Create plotly figures
fig1 <- plot_ly(x = x, y = y) %>%
add_surface(z = outer(x, y, f1), colorscale = 'Blues', name = 'x/y-1')
fig2 <- plot_ly(x = x, y = y) %>%
add_surface(z = outer(x, y, f2), colorscale = 'Reds', name = 'log(x)-log(y)')
fig3 <- plot_ly(x = x, y = y) %>%
add_surface(z = outer(x, y, f3), colorscale = 'Greens', name = 'x=y')
# Set axis labels and layout
fig1 <- fig1 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'z')))
fig2 <- fig2 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'z')))
fig3 <- fig3 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'z')))
# Arrange plots in a grid
subplot(fig1, fig2, fig3, nrows = 1)
## Warning: 'layout' objects don't have these attributes: 'NA'
## Valid attributes include:
## '_deprecated', 'activeshape', 'annotations', 'autosize', 'autotypenumbers', 'calendar', 'clickmode', 'coloraxis', 'colorscale', 'colorway', 'computed', 'datarevision', 'dragmode', 'editrevision', 'editType', 'font', 'geo', 'grid', 'height', 'hidesources', 'hoverdistance', 'hoverlabel', 'hovermode', 'images', 'legend', 'mapbox', 'margin', 'meta', 'metasrc', 'modebar', 'newshape', 'paper_bgcolor', 'plot_bgcolor', 'polar', 'scene', 'selectdirection', 'selectionrevision', 'separators', 'shapes', 'showlegend', 'sliders', 'smith', 'spikedistance', 'template', 'ternary', 'title', 'transition', 'uirevision', 'uniformtext', 'updatemenus', 'width', 'xaxis', 'yaxis', 'barmode', 'bargap', 'mapType'
# Display plot
fig1
library(plotly)
# Define x, y ranges
x <- seq(0.01, 10, length.out = 100)
y <- seq(0.01, 10, length.out = 100)
# Define functions
f1 <- function(x, y) x/(y+1) - log(x/y)
f2 <- function(x, y) x-y
# Create plotly figures
fig1 <- plot_ly(x = x, y = y) %>%
add_surface(z = outer(x, y, f1), colorscale = 'Blues', name = 'diff')
fig2 <- plot_ly(x = x, y = y) %>%
add_surface(z = outer(x, y, f2), colorscale = 'Greens', name = 'x=y')
# Set axis labels and layout
fig1 <- fig1 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'z')))
fig2 <- fig2 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'z')))
# Arrange plots in a grid
subplot(fig1, fig2, nrows = 1)
## Warning: 'layout' objects don't have these attributes: 'NA'
## Valid attributes include:
## '_deprecated', 'activeshape', 'annotations', 'autosize', 'autotypenumbers', 'calendar', 'clickmode', 'coloraxis', 'colorscale', 'colorway', 'computed', 'datarevision', 'dragmode', 'editrevision', 'editType', 'font', 'geo', 'grid', 'height', 'hidesources', 'hoverdistance', 'hoverlabel', 'hovermode', 'images', 'legend', 'mapbox', 'margin', 'meta', 'metasrc', 'modebar', 'newshape', 'paper_bgcolor', 'plot_bgcolor', 'polar', 'scene', 'selectdirection', 'selectionrevision', 'separators', 'shapes', 'showlegend', 'sliders', 'smith', 'spikedistance', 'template', 'ternary', 'title', 'transition', 'uirevision', 'uniformtext', 'updatemenus', 'width', 'xaxis', 'yaxis', 'barmode', 'bargap', 'mapType'
# Display plot
fig1
8.4 Pairs Trading
Implement a pairs trading strategy of your own, using daily stock data from Yahoo Finance. More specifically, find two ETF’s that are highly correlated (e.g. both track the same index or industry) and download 10 years of historical data. Use the first 5 years to select the parameters of your strategy, and the last 5 years to test it. Plot the log-difference of the series over the entire period, as well as the cumulative profit of the strategy
We are trading “SPY” against “IVV”; both are large ETFs tracking the S&P500 index. The plot of the log-ratio of their prices is shown below:
library(zoo)
library(tseries)
ETF.names=c('SPY','IVV')
N.ETF=length(ETF.names)
S=list()
for(i in 1:N.ETF){
print(paste(i, ' - ', ETF.names[i]))
S[[i]]=get.hist.quote(ETF.names[i], start='2006-01-03', end='2016-12-30', quote='AdjClose')
}
## [1] "1 - SPY"
## time series ends 2016-12-29
## [1] "2 - IVV"
## time series ends 2016-12-29
names(S)=ETF.names
n=length(S[[1]]); n.half=n%/%2
plot(log(S[[1]]/S[[2]]), type='l'); abline(v=index(S[[1]])[n.half], lty=2, lwd=2)
Note that the log-ratio exhibits highly mean reverting behavior in the training period (left half), especially around the Fall 2008 financial crisis. In the testing period (right half), there are less intense fluctuations, and an adjustment (jump) at the beginning of 2016. The pairs trading strategy employed involves opening a position at 2 standard deviations, closing a position at the mean, and bailing-out a position at 4 standard deviations (where the standard deviation is calculated based on the first half of data). The resulting strategy and cumulative gain (per $1 position in either stock) are shown below:
S.1.train=coredata(S[[1]][1:n.half])
S.2.train=coredata(S[[2]][1:n.half])
S.1.test=coredata(S[[1]][(n.half+1):n])
S.2.test=coredata(S[[2]][(n.half+1):n])
profit.PP=rep(0,length(S.1.test))
logRatio.train=log(S.1.train/S.2.train)
logRatio.test=log(S.1.test/S.2.test)
std.logRatio.test=(logRatio.test-mean(logRatio.train))/sd(logRatio.train)
open.lim=2
bail.lim=4
trade.ind=0 # trade indicator: 0 means no position, -1 means buy 1st & sell 2nd stock, 1 means sell 1st & buy 2nd stock
bail.ind=0 # bail-out indicator
t.open=t.close=NULL
for(j in 1:length(logRatio.test)){
# check for closing an open position
if( (trade.ind*std.logRatio.test[j]<0 | abs(std.logRatio.test[j])>bail.lim) & bail.ind==0 ){
trade.ind=0
profit.PP[j]=profit.PP[j]+n1*S.1.test[j]+n2*S.2.test[j]
t.close=c(t.close,j)
if( abs(std.logRatio.test[j])>bail.lim ){
bail.ind=1 # bail out position, & stop trading
}
}
# check for opening a position
if( (trade.ind==0) & (abs(std.logRatio.test[j])>open.lim) & (bail.ind==0) ){
trade.ind=sign(std.logRatio.test[j])
n1=-trade.ind/S.1.test[j]
n2=trade.ind/S.2.test[j]
t.open=c(t.open,j)
}
}
par(mfrow=c(2,1))
plot(cumsum(profit.PP), type='l', ylab="cumulative profit per $1 invested", lwd=2)
plot(std.logRatio.test, type='l', lwd=2, ylim=c(-4,4))
abline(h=c(-2,2), col=3); abline(h=c(-4,0,4), col=2)
abline(v=t.open, col=3, lty=2); abline(v=t.close, col=2, lty=2)
9. Monte Carlo Simulation
9.2
Consider the geometric Brownian motion \(S_t=50 \exp \left\{.03 t+.2 W_t\right\}\); fix \(n=100\) and let \(t_i=i /(n+1), \quad i=0, \ldots, n\). Write \(\mathrm{R}\) code that simulates and plots a conditional path of the process, sampled at times \(\left\{t_i\right\}_{i=0}^{n+1}\), given that \(W_1=1\).
Use the result of the previous question combined with the Cholesky decomposition for generating the correlated values of \(W\) at \(\left(t_1, \ldots, t_n\right)\).
set.seed(1234567890)
S0=50; mu=.03; sigma=.2; n=100; t.i=(1:n)/(n+1)
X1=mu+sigma*1; S1=S0*exp(X1)
MU.X=t.i*X1 # conditional mean
mat1=matrix(t.i,n,n)
mat2=matrix(t.i,n,n,byrow=T)
SIGMA.X= (pmin(mat1,mat2) - t.i %x% t(t.i))*sigma^2 # conditional covariance matrix
L=chol(SIGMA.X) # Cholesky decomposition of SIGMA.X
Z=rnorm(n)
X=MU.X + t(L)%*%Z
S=S0*exp(X)
plot(c(S0,S,S1), type='l')
9.3 Binary Option
Payoff: \(S_T \mathbf{1}_{\left\{S_T \geq K\right\}}=\left\{\begin{array}{ll}S_T, & S_T \geq K \\ 0, & S_T<K\end{array}\right.\) at maturity \(T\).
This is a European asset-or-nothing binary call, and has a analytical price in the Black-Scholes model (i.e. for constant interest rate and an underlying asset following Geometric Brownian motion).
Consider a model where the interest rate is \(r=.05\) and the asset has initial price \(S(0)=100\) and volatility \(\sigma=.4\). Find the price of the assetor-nothing call with strike \(K=100\) and maturity \(T=.5\) using Monte Carlo simulation with \(n=10,000\).
Report your point estimate and the approximate \(95 \%\) confidence interval for the price, and compare it to the exact price.
set.seed(1234567890)
S0=100; r=0.05; sigma=.4; T= 0.5; K = 100
# exact price
d1 = ( log(S0/K) + ( r + sigma^2 / 2 ) * T ) / ( sigma * sqrt(T) )
C.exact = S0 * pnorm( d1 )
print( paste("Exact binary call price =", C.exact) )
## [1] "Exact binary call price = 59.0880178044312"
# MC simulation
n=10000
Z = rnorm( n )
ST = S0 * exp( ( r - sigma^2 / 2 ) * T + sigma * sqrt(T) * Z )
payoff_3 = exp( -r * T ) * ifelse( ST>K, ST, 0 )
mean_3 = mean(payoff_3)
se_3 = sd(payoff_3) / sqrt(n)
CI_3 = mean_3 + c(-1,1) * qnorm( .975) * se_3
print( paste("MC binary call price =", mean_3) )
## [1] "MC binary call price = 58.2251962242492"
print( paste( c("MC binary call 95% CI =", CI_3), collapse = " " ) )
## [1] "MC binary call 95% CI = 56.984595876404 59.4657965720943"
9.4 Rainbow Option
A derivative exposed to two or more underlying assets.
Consider a European rainbow option with payoff given by \(\left(\max _i\left\{S_T^{(i)}\right\}-\frac{1}{d} \sum_{i=1}^d S_T^{(i)}\right)\), i.e. the payoff is the maximum final price minus the average final price of all \(d\) assets (note that the maximum/average is over assets, not time). Estimate the price of this option using Monte Carlo simulation with \(n=10,000\) d-dimensional paths, and provide a \(95 \%\) confidence interval with your answer. Assume that \(d=5, K=100, T=1, r=0.03, S^{(i)}(0)=\) 100, \(\forall i=1, \ldots, d\), and that the assets follow multivariate Geometric Brownian Motion: \[ \begin{aligned} & d \mathbf{S}(t)=r \mathbf{S}(t) d t+\sigma \circ \mathbf{S}(t) \circ d \mathbf{W}(t) \\ & \text { for } \sigma=\left[\begin{array}{lllll} .1 & .2 & .3 & .4 & .5 \end{array}\right] \text { and } \operatorname{Corr}\left(W_i(1), W_j(1)\right)=.3, \quad \forall i \neq j \end{aligned} \] where \(\circ\) is the Hadamard or element-wise matrix product. (Note: there is no general analytical solution for rainbow options.)
set.seed(1234567890)
n=10000; d=5; S0=100; M=1; r=.01
Rho= matrix(.3,5,5); diag(Rho)=1 # correlation matrix
V=c(1:5/10) # volatilities
Sig=Rho*(V%*%t(V)) # covariance matrix
L=chol(Sig) # Cholesky factorization of covariance matrix
Drft=matrix(r-V^2/2, n, d, byrow=TRUE) # drift
Z=matrix( rnorm(n*d) , n, d)
ST=S0*exp( Drft*M+ Z%*%L*sqrt(M) )
maxT=apply(ST,1,max)
avgT=apply(ST,1,mean)
payoff_4=exp(-r*M)*(maxT-avgT)
mean_4=mean(payoff_4)
se_4=sd(payoff_4)/sqrt(n)
CI_4=mean_4+c(-1,1)*qnorm(.975)*se_4
print( paste("rainbow option price =", mean_4) )
## [1] "rainbow option price = 36.2473636425247"
print( paste( c("rainbow option 95% CI =", CI_4), collapse = " " ) )
## [1] "rainbow option 95% CI = 35.7261581492228 36.7685691358267"
10. Simulation: Pricing Exotic Derivatives
10.2 Extrema of BM
A European lookback option is a path-dependent option whose payoff at maturity depends on the maximum/minimum price of the underlying asset before maturity. There are two types of lookback options, namely fixed strike & floating strike, with payoffs:
| Call | Put | |
|---|---|---|
| fixed strike | \((M_T - K)_+\) | \((K-m_T)_+\) |
| float strike | \((S_T - m_T)_+\) | \((M_T - S_T)_+\) |
where \(M_T = \underset{0\leq t \leq T}\max\{S_t\}\) and \(m_T = \underset{0\leq t \leq T}\min\{S_t\}\)
(a) floating strike
Perform Monte Carlo simulation for estimating the price of a floating strike European lookback call option. Assume the underlying asset price follows Geometric BM with \(S_0=90, \quad T=1, r=2 \%, \sigma=20 \%\), and use unbiased estimation (i.e. simulate from the exact distribution of the minimum). Use \(n=100,000\) samples and create a \(95 \%\) confidence interval. Compare your result with the exact price of the option, which is 14.26674
set.seed(12345)
n=1000000 # number of paths
r=.02 # risk-free rate
v=.2 # volatility
M=1 # maturity
S0=90 # starting asset price
# use reflection trick for simulating minimum:
# minimum of arithmetic BM is minus the maximum of process with opposite drift
Z=rnorm(n)
XT.reflected=-(r-v^2/2)*M + v*sqrt(M)*Z # use negative drift for reflected path
MT.reflected=(XT.reflected + sqrt( XT.reflected^2 -2*v^2*M*log(runif(n))))/2 # reflected max
ST=S0*exp(-XT.reflected) # final value
mT=S0*exp(-MT.reflected) # minimum
payoff_2.a = exp(-r*M) * (ST-mT)
mean_2.a=mean(payoff_2.a) # MC estimate
se_2.a=sd(payoff_2.a)/sqrt(n) # st. dev
CI_2.a=mean_2.a+c(-1,1)*qnorm(.975)*se_2.a # 95% CI
print( paste("MC Float Strike Lookback call price =", mean_2.a) )
## [1] "MC Float Strike Lookback call price = 14.2733852982329"
print( paste( c("MC Float Strike Lookback call 95% CI =", CI_2.a), collapse = " " ) )
## [1] "MC Float Strike Lookback call 95% CI = 14.2487087208825 14.2980618755833"
print( paste("Exact Float Strike Lookback call price =", 14.2667437) )
## [1] "Exact Float Strike Lookback call price = 14.2667437"
(b) fixed strike
Perform Monte Carlo simulation for estimating the price of a fixed strike European lookback call option. Assume the underlying asset price follows Geometric BM with \(S_0=90, \quad K=90, \quad T=1, r=2 \%, \sigma=20 \%\), and use unbiased estimation (i.e. simulate from the exact distribution of the maximum). Use \(n=100,000\) samples and create a \(95 \%\) confidence interval. Compare you result with the exact price of the option, which is 16.04886
K=90 # fixed strike
Z=rnorm(n)
XT=(r-v^2/2)*M + v*sqrt(M)*Z
MT=S0 * exp ( (XT+sqrt(XT^2 -2*v^2*M*log(runif(n))))/2 )
payoff_2.b = exp(-r*M)*pmax(MT-K,0)
mean_2.b=mean(payoff_2.b) # MC estimate
se_2.b=sd(payoff_2.b)/sqrt(n) # st. dev
CI_2.b=mean_2.b+c(-1,1)*qnorm(.975)*se_2.b # 95% CI
print( paste("MC Fixed Strike Lookback call price =", mean_2.b) )
## [1] "MC Fixed Strike Lookback call price = 16.0314639209374"
print( paste( c( "MC Fixed Strike Lookback call 95% CI =", CI_2.b), collapse = " " ) )
## [1] "MC Fixed Strike Lookback call 95% CI = 16.0052209559405 16.0577068859343"
print( paste("Exact Fixed Strike Lookback call price =", 16.0488631) )
## [1] "Exact Fixed Strike Lookback call price = 16.0488631"
(c) Euler discretization
Repeat the previous part assuming the risk-neutral asset price dynamics: \[ d S_t=\mu\left(t, S_t\right) d t+\sigma\left(t, S_t\right) d W_t=\left(r S_t\right) d t+\left(\sigma \cos \left(e^t\right) \sqrt{S_t}\right) d W_t \] Note that the price process is not a Geometric Brownian Motion anymore. Use Euler discretization with \(m=50\) steps and simulate \(n=10,000\) paths. In order to approximate the maximum of each path, simulate the maximum of each step using the result for the extrema of an arithmetic BM: \[ M_{t_i} \mid S_{t_{i-1}}, S_{t_i}=S_{t_{i-1}}+\frac{\Delta S_{t_i}+\sqrt{\Delta S_{t_i}^2-2 \cdot \Delta t \cdot \sigma^2\left(t_{i-1}, S_{t_{i-1}}\right) \cdot \log (U)}}{2} \] Where \[ M_{t_i}=\max _{t_{i-1} \leqslant t \leqslant t_i}\left\{S_t\right\}, \Delta S_{t_i}=S_{t_i}-S_{t_{i-1}}, \Delta t=t_i-t_{i-1} \& U \sim \text { Uniform }(0,1) \] The maximum of the entire path will be the maximum over all steps in the path, i.e. \(M_T=\max _{0 \leqslant t \leqslant T}\left\{S_t\right\}=\max _{i=0, \ldots, m}\left\{M_{t_i}\right\}\). Include a \(95 \%\) confidence interval with your answer.
K=90
n=10000; m=50; Dt=M/m
ST=matrix(S0,n,m+1); MT=matrix(0,n,m)
for(i in 1:m){
ST.prev=ST[,i]
step.mu=r*ST.prev*Dt;
step.sig=v*cos(exp(i*Dt))*sqrt(ST.prev)*sqrt(Dt)
DST=step.mu+step.sig*rnorm(n)
ST[,i+1]=ST[,i]+DST
MT[,i]=ST[,i]+(DST+sqrt(DST^2-2*step.sig^2*log(runif(n))))/2
}
MT.all=apply(MT,1,max)
ST.final=ST[,m+1]
payoff_2.c = exp(-r*M)*pmax(MT.all-K,0)
mean_2.c=mean(payoff_2.c) # MC estimate
se_2.c=sd(payoff_2.c)/sqrt(n) # st. dev
CI_2.c=mean_2.c+c(-1,1)*qnorm(.975)*se_2.c # 95% CI
print( paste("MC Fixed Strike Lookback call price =", mean_2.c) )
## [1] "MC Fixed Strike Lookback call price = 2.18656942157468"
print( paste( c( "MC Fixed Strike Lookback call 95% CI =", CI_2.c), collapse = " " ) )
## [1] "MC Fixed Strike Lookback call 95% CI = 2.17236758314114 2.20077126000821"
10.3 European rainbow option
Its payoff is given by \(\left(\max _i\left\{S_T^{(i)}\right\}-\frac{1}{d} \sum_{i=1}^d S_T^{(i)}\right)\), i.e. the payoff is the maximum final price minus the average final price of all \(d\) assets (note that the maximum/average is over assets, not time).
Estimate the price of this option using Monte Carlo simulation with \(n=10,000 \mathrm{~d}\)-dimensional paths, and provide a \(95 \%\) confidence interval with your answer.
Assume that \(d=5, K=100, T=1, r=0.03, S^{(i)}(0)=100, \forall i=1, \ldots, d\), and that the assets follow the multivariate SDE: \[ d \mathbf{S}(t)=r \circ \mathbf{S}(t) d t+\sigma \circ \breve{\mathbf{S}}(t) \circ d \mathbf{W}(t) \] for \(\sigma=\left[\begin{array}{lllll}.1 & .2 & .3 & .4 & .5\end{array}\right], \operatorname{Corr}\left(W_i(1), W_j(1)\right)=.3, \quad \forall i \neq j\)
and \(\breve{\mathbf{S}}(t)=\left[\begin{array}{lll}S^{(5)}(t) & \cdots & S^{(1)}(t)\end{array}\right]^T\) is the reverse of \(\mathbf{S}(t)\)
This process does not follow a multivariate Geometric BM, so use Euler discretization with \(m=25\) steps to approximate the final prices of the assets.
d = 5 #dimension
Rho= matrix(.3,d,d); diag(Rho)=1 # correlation matrix
V=c(1:5/10) # volatilities
Sig=Rho*(V%*%t(V)) # covariance matrix
L=chol(Sig) # Cholesky factorization of covariance matrix
Drft=matrix(r-V^2/2, n, d, byrow=TRUE) # drift
m=25; Dt=M/m
ST=matrix(S0, n, d)
for(i in 1:m){
Z=matrix(rnorm(n*d),n,d)
ST=ST+(r*ST*Dt + Z%*%L*sqrt(Dt)*ST[,d:1] )
ST[ST<0]=0 # set possible negative prices to 0
}
maxT=apply(ST,1,max)
avgT=apply(ST,1,mean)
payoff_3=exp(-r*M)*(maxT-avgT)
mean_3=mean(payoff_3)
se_3=sd(payoff_3)/sqrt(n)
CI_3=mean_3+c(-1,1)*qnorm(.975)*se_3
print( paste("rainbow option price =", mean_3) )
## [1] "rainbow option price = 30.56361486597"
print( paste( c("rainbow option 95% CI =", CI_3), collapse = " " ) )
## [1] "rainbow option 95% CI = 30.2443640059131 30.882865726027"
11. Simulation: Variance Reduction Techniques
set.seed(12345)
# Function to compute Black-Scholes prices of European call & put options
BSprice=function(S, X, r, M, v){
# arguments: Asset price (S), Strike (X), risk-free rate (r)
# Maturity (M), volatility (v)
d1=(log(S/X)+(r+0.5*v^2)*M)/(v*sqrt(M))
d2=d1-v*sqrt(M)
call.pr = S*pnorm(d1) - X*exp(-r*M)*pnorm(d2)
put.pr = X*exp(-r*M) * pnorm(-d2) - S*pnorm(-d1)
return(cbind(call.pr,put.pr))
}
11.2 Gap Option
A European gap option has a strike price K1 and a trigger price K2. The trigger price determines whether or not the gap option is exercised, while the strike price determines the amount of the payoff. For example, the European gap call option has payoff given by:
\[ \begin{cases} S_T - K_1 & S_T \geq K_2\\ 0 & S_T < K_2 \end{cases} \]
Obviously, if the strike price is equal to the trigger price then the gap option is an ordinary option. Note that gap options can result in negative payoff at expiration (so they should not really be called “options”). Moreover, gap options are different from barrier options, in that for barrier options the barrier can be reached prior to expiration, whereas for gap options only the final stock price is compared to the trigger.
(a) GBM price
Estimate the price of a European gap call option with \(S_0=55, K_1=50, K_2=60, T=1 y r, r=2 \%, \sigma=20 \%\) and assuming the price follows GBM. Use \(n=100,000\) paths and calculate the standard error of your result. Compare that to the exact price of 6.129964.
n=100000; S0=55; K1=50; K2=60; M=1; r=.02; v=.2
Z=rnorm(n)
ST=S0*exp( (r-v^2/2)*M + v*sqrt(M)*Z)
payoff_2.a=exp(-r*M)*(ST-K1)*(ST>K2)
(mean_2.a=mean(payoff_2.a))
## [1] 6.121
(se_2.a=sd(payoff_2.a)/sqrt(n))
## [1] 0.03081
(b) antithetic variables
Repeat part (a) using antithetic variables, and report the new price estimate and standard error.
Z=rnorm(n/2); Z=cbind(Z,-Z)
ST=S0*exp( (r-v^2/2)*M + v*sqrt(M)*Z)
payoff_2.b=rowMeans(exp(-r*M)*(ST-K1)*(ST>K2))
(mean_2.b=mean(payoff_2.b))
## [1] 6.121
(se_2.b=sd(payoff_2.b)/sqrt(n))
## [1] 0.01694
(c) stratification
Repeat part (a) using stratification with m = 5, 10, 20 equiprobable strata, and report the new price estimate and standard error
m.all=c(5,10,20)
mean_2.c=se_2.c=rep(0,length(m.all))
# iterate over 3 strata values
for(i in 1:length(m.all)){
m=m.all[i]
nj=n/m # num paths / num strata
Z=matrix(0,nj,m)
# iterate over each of m strata
for(j in 1:m){
temp=runif(nj,(j-1)/m,j/m)
Z[,j]=qnorm(temp)
}
ST=S0*exp( (r-v^2/2)*M + v*sqrt(M)*Z )
payoff_2.c=exp(-r*M)*(ST-K1)*(ST>K2)
mean_2.c[i]=mean( colMeans(payoff_2.c) )
se_2.c[i]=sqrt( sum( apply(payoff_2.c,2,var)/nj ) )/m
}
cbind(mean_2.c, se_2.c)
## mean_2.c se_2.c
## [1,] 6.124 0.013190
## [2,] 6.129 0.008517
## [3,] 6.131 0.005919
(d) comparison
Would you expect the variance reduction methods in part (b) and (c) to work better for in-the-money, at-the-money, or out-of-the-money options? Justify your answer
Antithetic variables: works best when original and antithetic variates are negatively related
Stratification: works best for highly variable payoffs (if target RV changes over its domain)
11.3 Forward start option
A forward start option is an option that starts at a specified future date with an expiration date set even further in the future. For example, consider a forward start call option with start date in \(1 \mathrm{yr}\) and expiration date in \(2 \mathrm{yrs}\) : at \(t=0\) you pay the option premium, and at \(t=1\) you receive a European call option with expiration date \(t=2\) and strike price \(K=S_1\) (this is called an at-the-money forward start call, because the strike price is set to be the asset price at \(t=1\), which is unknown at \(t=0\) ).
(a) GBM price
Estimate the price of an at-the-money forward start European call option with \(T_{\exp }=2 y r, r=2 \%, \sigma=20 \%\) and assuming the price follows GBM. Use \(n=10,000\) paths and calculate the standard error of your result. You should use the BlackScholes formula to explicitly calculate the value of the option at \(T_{\text {start }}=1\), i.e. you should generate \(S_1\) and then calculate its “payoff” as the Black-Scholes price of a European at-the-money call option with 1 year to maturity (use the BSprice() function provided in EX9.R). Compare that to the exact price of 6.241225.
n=10000; S0=70; Tstart=1; Texp=2; r=.02; v=.2
Z=rnorm(n)
STstart=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z )
payoff_3.a = exp(-r*Tstart) * BSprice(STstart, STstart, r, Texp-Tstart, v)[,1]
(mean_3.a=mean(payoff_3.a)) # MC estimate
## [1] 6.246
(se_3.a=sd(payoff_3.a)/sqrt(n))
## [1] 0.01266
(b) antithetic variables
Repeat part (a) using antithetic variables, and report the new price estimate and standard error.
Z=rnorm(n/2); Z.a=-Z
STstart=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z )
STstart.a=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z.a ) # antithetic path
payoff=exp(-r*Tstart) * BSprice(STstart, STstart, r, Texp-Tstart, v)[,1]
payoff.a=exp(-r*Tstart) * BSprice(STstart.a, STstart.a, r, Texp-Tstart, v)[,1] # antithetic path payoff
payoff_3.b=(payoff+payoff.a)/2
(mean_3.b=mean(payoff_3.b)) # MC estimate
## [1] 6.244
(se_3.b=sd(payoff_3.b)/sqrt(length(payoff_3.b))) # std. error
## [1] 0.002629
(c) stratification
Repeat part (a) using stratification with equiprobable strata, and report the new price estimate and standard error.
m=10; nj=n/m
Z=STstart=payoff_3.c=matrix(0,nj,m)
for(j in 1:m){
temp=runif(nj,(j-1)/m,j/m)
Z[,j]=qnorm(temp)
STstart[,j]=S0*exp( (r-v^2/2)*Tstart + v*sqrt(Tstart)*Z[,j] )
payoff_3.c[,j]=exp(-r*Tstart) * BSprice(STstart[,j], STstart[,j], r, Texp-Tstart, v)[,1]
}
(mean_3.c=mean( colMeans(payoff_3.c) ) )
## [1] 6.239
(se_3.c=sqrt( sum( apply(payoff_3.c,2,var)/nj ) )/m )
## [1] 0.002762
(d) proof
(paper & pencil part) Show that the price of the at-the-money forward call option is actually a simple linear function of \(S_0\)
(Hint: use risk-neutral pricing and write the payoff of the option at \(T_{\text {start }}\) as the Black-Scholes price of an at-the-money call with maturity at \(T>T_{\text {start }}\) )
Dt=Texp-Tstart
d1=(r/v+v/2)*sqrt(Dt)
d2=(r/v-v/2)*sqrt(Dt)
(exact_price=S0*( pnorm(d1) - exp(-r*Dt)*pnorm(d2) ))
## [1] 6.241
The Black-Scholes price for a European call is \[ C=S_0 \Phi\left(d_1\right)-e^{-r T} K \Phi\left(d_2\right) \] where \(d_1=\frac{\ln \left(S_0 / K\right)+\left(r+\frac{1}{2} \sigma^2\right) T}{\sigma \sqrt{T}} ; d_2=d_1-\sigma \sqrt{T}\).
For an at-the-money call, in particular, we get \(C=S_0\left[\Phi\left(d_1\right)-e^{-r T} \Phi\left(d_2\right)\right]=\) \(S_0\left[\Phi\left(\frac{\left(r+\frac{1}{2} \sigma^2\right) T}{\sigma \sqrt{T}}\right)-e^{-r T} \Phi\left(\frac{\left(r-\frac{1}{2} \sigma^2\right) T}{\sigma \sqrt{T}}\right)\right]\)
By risk-neutral pricing, the price of the at-the-money forward call option becomes: \[ \begin{aligned} C_{F W D} & =\mathbb{E}\left[e^{-r T_{\exp }} \text { payoff }\right] \\ & =\mathbb{E}\left\{\mathbb{E}\left[e^{-r T_{\exp }}\left(\text { at the money Euro call at } T_{\text {start }}\right) \mid S_{T_{\text {start }}}\right]\right\} \\ & =\mathbb{E}\left[e^{-r T_{\text {start }}} C\left(K=S_{T_{\text {start }}}, T=T_{\exp }-T_{\text {start }}\right)\right]\\ & =\alpha \mathbb{E}\left[e^{-r T_{\text {start }}} S_{T_{\text {start }}}\right]=\alpha S_0 \end{aligned} \] where \[ \alpha = \Phi\left(\frac{\left(r+\frac{1}{2} \sigma^2\right)\left(T_{\exp }-T_{\text {start }}\right)}{\sigma \sqrt{T_{\exp }-T_{\text {start }}}}\right)-e^{-r\left(T_{\exp }-T_{\text {start }}\right)} \Phi\left(\frac{\left(r-\frac{1}{2} \sigma^2\right)\left(T_{\exp }-T_{\text {start }}\right)}{\sigma \sqrt{T_{\exp }-T_{\text {start }}}}\right) \] In the last line we used the fact that the discounted asset price is a martingale under the risk-neutral measure.
11.4 alpha quantile option
An \(\alpha\)-quantile option is an option whose payoff is determined by the \(\alpha\) th quantile of the underlying price from time 0 until expiration. For example, the median call option (i.e., \(\alpha=50 \%)\) with strike \(K\) has payoff \(\left(\operatorname{median}\left\{S_t\right\}_{t=0}^T-K\right)_{+}\), where median \(\left\{S_t\right\}_{t=0}^T\) is the median asset price over \([0, T]\). Assume that the underlying asset follows GBM with \(S_0=60, r=5 \%, \sigma=25 \%\), and use simulation to estimate the price of median call with \(K=70 \& T=1\)
(a) path discretization
Perform simple Monte Carlo using path discretization with \(t_j=j \frac{T}{m}, j=0, \ldots, m\) and \(m=50\) in order to approximate the continuous median \(\left\{S_t\right\}_{t=0}^T\) by the discrete median \(\left\{S\left(t_j\right)\right\}_{j=0}^m\). Use \(n=10,000\) paths and report the standard error of the estimate.
S0=60; K=70; M=1; r=.05; v=.25
n=10000; m=50; Dt=M/m
Z=matrix(rnorm(n*m), n, m)
S=S0*exp(t(apply((r-v^2/2)*Dt+Z*v*sqrt(Dt),1,cumsum)))
medS=apply(S,1,median)
payoff_4.a=exp(-r*M)*pmax(medS-K,0)
(mean_4.a=mean(payoff_4.a))
## [1] 1.036
(se_4.a=sd(payoff_4.a)/sqrt(n))
## [1] 0.03346
(b) control variates
Repeat part (a) using the Black-Scholes price of the European call as a control variate, and report the new price estimate and standard error
EX=BSprice(S0, K, r, M, v)[1]
X=exp(-r*M)*pmax(S[,m]-K,0) # control variate values
Y=payoff_4.a
(rho_XY=cor(X,Y))
## [1] 0.7177
(b=sd(Y)/sd(X)*rho_XY)
## [1] 0.3005
(mean_4.b=mean_4.a - b*(mean(X)-EX))
## [1] 1.039
(se_4.b=se_4.a*sqrt((1-rho_XY^2)))
## [1] 0.0233
11.5 Asian option
An Asian option is an option whose payoff is determined by the average underlying price over some preset period of time, typically the time till expiration. For example, the Asian call option with strike \(K\) has payoff \((A(0, T)-K)_{+}\), where \(A(0, t)\) is the average asset price over \((0, T)\). There are two types of averaging: \[ \left\{\begin{array}{l} \text { arithmetic average : } A(0, T)=\frac{1}{T} \int_0^T S(t) d t \\ \text { geometric average : } A(0, T)=\exp \left(\frac{1}{T} \int_0^T \ln (S(t)) d t\right) \end{array}\right. \] If the underlying asset follows GBM, then the price of an Asian option with geometric averaging can be calculated explicitly. For Asian options with arithmetic averaging, however, we rely exclusively on numerical methods for pricing. Assume that the underlying asset follows GBM with \(S_0=70, r=5 \%, \sigma=25 \%\), and use simulation to estimate the price of an arithmetic Asian call with \(K=80 \& T=1\).
set.seed(12345)
S0=70; K=80; M=1; r=.05; v=.25
n=10000; m=50; Dt=M/m
(a) path discretization
Perform simple Monte Carlo using path discretization with \(m=50\) in order to approximate \[ A(0, T)=\frac{1}{T} \int_0^T S(t) d t \] by \(\tilde{A}(0, T)=\frac{1}{m} \sum_{j=1}^m S\left(t_i\right)\). Use \(n=10,000\) paths and report the standard error of the estimate.
Z=matrix(rnorm(n*m), n, m)
S=S0*exp(t(apply((r-v^2/2)*Dt+Z*v*sqrt(Dt),1,cumsum)))
A=rowMeans(S)
payoff_5.a=exp(-r*M)*pmax(A-K,0)
(mean_5.a=mean(payoff_5.a))
## [1] 1.493
(se_5.a=sd(payoff_5.a)/sqrt(n))
## [1] 0.04219
(b) control variate
Repeat part (a) using the price of the geometric Asian call as a control variate, and report the new price estimate and standard error. For implementing the control variate, the exact price of the geometric Asian option is 1.271969.
EX=1.271969 # Control variate - exact Geom. Asian call price
A.geom=exp(rowMeans(log(S)))
X=exp(-r*M)*pmax(A.geom-K,0) # control variate values
Y=payoff_5.a
(rho_XY=cor(X,Y))
## [1] 0.9987
(b=sd(Y)/sd(X)*rho_XY)
## [1] 1.081
(mean_5.b=mean_5.a - b*(mean(X)-EX))
## [1] 1.405
(se_5.b=se_5.a*sqrt((1-rho_XY^2)))
## [1] 0.002139
(c) importance sampling
Repeat part (a) using importance sampling. Change the probability measure so that \(E\left[S_T\right]=K\), and experiment with at least two more measure changes of your choice. Report the new price estimates and standard errors.
MEANS=c(80,90,100)
mean_5.c=se_5.c=rep(0,length(MEANS))
for(i in 1:length(MEANS)){
new.r=log(MEANS[i]/S0)/M
ZZ=(new.r-v^2/2)*Dt+Z*v*sqrt(Dt)
S=S0*exp(t(apply(ZZ,1,cumsum)))
A=rowMeans(S)
log.ratio=dnorm(ZZ,(r-v^2/2)*Dt,v*sqrt(Dt),log=TRUE)-dnorm(ZZ,(new.r-v^2/2)*Dt,v*sqrt(Dt),log=TRUE)
prod.ratio=exp( rowSums(log.ratio) )
payoff=exp(-r*M)*pmax(A-K,0)*prod.ratio
mean_5.c[i]=mean(payoff)
se_5.c[i]=sd(payoff)/sqrt(n)
}
cbind(mean_5.c, se_5.c)
## mean_5.c se_5.c
## [1,] 1.480 0.02973
## [2,] 1.482 0.02138
## [3,] 1.494 0.01903
12. Optimization in Finance
12.1
Consider the example of trading an asset with transaction costs, assuming you have perfect knowledge of its price. Adapt the code from lecture 12 to find the optimal strategy and its cumulative P/L, assuming you can also short the asset, i.e. you have to consider three states: long, neutral, and short. Note that when you go from long to short position, or vice-versa, you have to pay twice the transaction cost because you sell/buy two units of asset. Test your code on Shopify (SHOP.TO) closing daily prices from Jan 1 to Mar 31, 2023 with a transaction cost of $0.1/share, and plot the strategy.
Using Dynamic Programming, find the optimal trading strategy under transaction costs (\(tc\)) and a given price evolution.
library(tseries)
St = get.hist.quote( "SHOP.TO", start = "2023-01-01", end = "2023-03-31", quote = "Close" )
## time series starts 2023-01-03
## time series ends 2023-03-30
plot(St, lwd=2, type = "o", pch=16, main = "Shopify Close Price")
tc = 1 # cost per transaction
n = length(St)
# optimal value function
n.val = rep( 0, n ) # neutral (n)
l.val = rep( 0, n ) # long (l)
s.val = rep( 0, n ) # short (s)
l.val[n] = as.numeric(St[n])-tc
s.val[n] = -as.numeric(St[n])-tc
# optimal strategy
n.str = rep( 0, n )
l.str = rep( 1, n )
s.str = rep( -1, n )
# Backward Induction/Dynamic Programming
# We have 3 options/transitions to consider for each state:
for( i in (n-1):1 ){
# entering at n, compare: n2n, n2l, n2s
n2n = n.val[i+1]
n2l = l.val[i+1] - as.numeric(St[i]) - tc
n2s = s.val[i+1] + as.numeric(St[i]) - tc
if( n2n >= max( n2l, n2s) ){ # action: n -> n
n.val[i] = n2n
n.str[i] = 0
}else if( n2l >= max( n2n, n2s) ){ # action: n -> l
n.val[i] = n2l
n.str[i] = 1
n.str[(i+1):n] = l.str[(i+1):n]
}else{ # action: n -> s
n.val[i] = n2s
n.str[i] = -1
n.str[(i+1):n] = s.str[(i+1):n]
}
# entering at l, compare: l2n, l2l, l2s
l2n = n.val[i+1] + as.numeric(St[i]) - tc
l2l = l.val[i+1]
l2s = s.val[i+1] + 2*( as.numeric(St[i]) - tc )
if( l2n >= max( l2l, l2s ) ){ # action: l -> n
l.val[i] = l2n
l.str[i] = 0
l.str[(i+1):n] = n.str[(i+1):n]
}else if( l2l >= max( l2n, l2s ) ){ # action: l -> l
l.val[i] = l2l
}else{
l.val[i] = l2s # action: l -> s
l.str[i] = -1
l.str[(i+1):n] = s.str[(i+1):n]
}
# entering at s, compare: s2n, s2l, s2s
s2n = n.val[i+1] - as.numeric(St[i]) - tc
s2l = l.val[i+1] + 2*(-as.numeric(St[i]) - tc)
s2s = s.val[i+1]
if( s2n >= max( s2l, s2s ) ){ # action: s -> n
s.val[i] = s2n
s.str[i] = 0
s.str[(i+1):n] = n.str[(i+1):n]
}else if( s2l >= max( s2n, s2s ) ){ # action: s -> l
s.val[i] = s2l
s.str[i] = 1
s.str[(i+1):n] = n.str[(i+1):n]
}else{
s.val[i] = s2s # action: s -> s
s.str[i] = -1
}
}
# calculate cumulative optimal value
run.val = rep(0, n)
run.val[1] = - n.str[1] * as.numeric(St[1]) - abs(n.str[1]) * tc
for(i in 2:n){
if( n.str[i] == n.str[i-1] ){
run.val[i] = run.val[i-1]
}else if( n.str[i] > n.str[i-1] ){
run.val[i] = run.val[i-1] + (n.str[i] - n.str[i-1]) * (-as.numeric(St[i]) - tc)
}else{
run.val[i] = run.val[i-1] + (n.str[i-1] - n.str[i]) * (as.numeric(St[i]) - tc)
}
}
run.val[n]
## [1] 58.52
n.val[1]
## [1] 58.52
# create time series
n.val = zoo::zoo( n.val, order.by = time(St))
n.str = zoo::zoo( n.str, order.by = time(St))
run.val = zoo::zoo( run.val, order.by = time(St))
# plots
tmp = 1 + n.str * 2
tmp = replace( tmp, n.str == -1, 2 )
plot( St, lwd=2, main = "Strategy")
points( St, col = tmp, pch = 16)
abline( v = time(St), lwd=.1, lty=2)
legend( "topleft", pch=16, col = c(2,1,3), c("short","neutral","long"))
plot( run.val, lwd=2 , main = "Value of strategy")
points( run.val, col = tmp, pch = 16)
abline( v = time(St), lwd=.1, lty=2)
legend( "topleft", pch=16, col = c(2,1,3), c("short","neutral","long"))
Finals
2021 Q3 Simulation & Betting Strategies
Consider the same setup as in question 1., i.e. a sequence of investment opportunities with i.i.d. Uniform \((-1,2)\) returns and an initial capital of \(V_0=1\).
(a) [5 points]
Use the optimize () function to find the Kelly fraction \(f^*\) to invest, i.e. value of \(f\) which maximizes the expression in Q1.(g).
E_log_wealth = function(x){
return( (2+1/x)*log(1+2*x) - (1/x-1)*log(1-x) - 1 )
}
f_ = optimize(E_log_wealth, c(0, 1), tol = 0.0001, maximum = T)$maximum
x = seq(0,1,.01)
plot(x, E_log_wealth(x), type ="l" )
points( f_, E_log_wealth(f_), pch = 16)
(b) [10 points]
Simulate 1,000 wealth paths consisting of 30 investments steps each using the Kelly fraction from the previous part (or, if you didn’t get it, just \(f=.7\) ). On the same axes, plot the mean simulated wealth, and the \(5 \%\) and \(95 \%\) quantiles of the simulated wealth at each step.
set.seed(123)
n=30; m=1000
R = matrix( runif(n*m, -1,2), n, m)
V = apply( 1+f_*R, 2, cumprod )
Vavg = apply( V, 1, mean )
Vlo = apply( V, 1, quantile, prob = .05 )
Vhi = apply( V, 1, quantile, prob = .95 )
plot( Vavg, type = "l", ylim = c(.001,25000) )
lines( Vlo, col=2 ); lines( Vhi, col = 2 )
# in log-scale
plot( Vavg, type = "l", log = "y", ylim = c(.001,25000) )
(c) [10 points]
In class, we mentioned that the Kelly criterion strategy reaches any wealth level faster than other strategies. You will test this claim with a simulation experiment. Fix a target wealth of 10,000 and simulate 1,000 paths with as many steps necessary to reach or exceed this level, under both the Kelly strategy \(\left(f^*\right)\) and the all-in strategy \((f=1)\). Report the average number of steps required to reach or exceed the target wealth under each strategy. Does it confirm the claim?
T_Kelly = T_all_in = integer(1000)
set.seed(123)
for( i in 1:1000 ){
# for Kelly
V = 1; j = 0
while( V < 10000 ){
V = V * (1 + f_ * runif(1,-1,2) )
j = j + 1
}
T_Kelly[i] = j
# for all-in
V = 1; j = 0
while( V < 10000 ){
V = V * (1 + runif(1,-1,2) )
j = j + 1
}
T_all_in[i] = j
}
mean(T_Kelly)
## [1] 54.31
mean(T_all_in)
## [1] 99.17
The Kelly criterion takes a shorter time, on average, to reach the goal
2021 Q4 Industry factor Models
We have seen how to fit macroeconomic factor models, by regressing returns \((\boldsymbol{R})\) on observed factors \((\boldsymbol{F})\) one asset at a time, to estimate their loadings \((\boldsymbol{\beta})\). There is a related class of models called fundamental factor models, which are also fit by regression, but taking a different approach. The idea is to look at asset characteristics (i.e. fundamentals) to determine the asset loadings \((\boldsymbol{\beta})\) based of their grouping with respect to common but unobserved factors. The simplest method is to group assets in non-overlapping industry categories, and assume all assets in a category are exposed to the same factor. In terms of the model, their loadings \((\beta\) ’s) are either 1 for their industry’s factor, or 0 for all other industry factors. The following model with 4 assets in 2 industries serves as an illustration: \[ \left[\begin{array}{l} R_{1, t} \\ R_{2, t} \\ R_{3, t} \\ R_{4, t} \end{array}\right]=\left[\begin{array}{ll} 1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 1 \end{array}\right]\left[\begin{array}{l} F_{1, t} \\ F_{2, t} \end{array}\right]+\left[\begin{array}{l} e_{1, t} \\ e_{2, t} \\ e_{3, t} \\ e_{4, t} \end{array}\right] \Leftrightarrow \boldsymbol{R}=\boldsymbol{\beta}^{\top} \boldsymbol{F}+\boldsymbol{e} \] In order to estimate the factor values \(\left(F_{1 / 2, t}\right)\), we use cross-sectional regression, i.e. we fix time \(t\) and look at returns as observations from the model \[ R_{i, t}=F_{1, t} I_{i, 1}+F_{2, t} I_{i, 2}+e_{i, t}, i=1, \ldots, n \] where \(I_{i, 1 / 2}\) are the indicator variables of asset \(i\) belonging to industry \(1 / 2\), serving as the explanatory variables, and the factors values \(F_{1 / 2, t}\) are the regression coefficients 1 . The following code loads 250 daily net returns for \(n=121\) assets (matrix \(\mathrm{R}\)) belonging in 2 industries (factor sector). The loop fits 250 regressions to estimate the factor values (matrix \(f\) ) and asset idiosyncratic errors (matrix e).
load("/Users/summer/Library/CloudStorage/GoogleDrive-summeryang720@gmail.com/My Drive/Obsidian/5B/STAD70/Finals/STAD70_W21_Final.RData")
b = model.matrix( object = ~ sectors - 1 ) # design matrix (betas)
f = matrix(0, nrow = nrow(R), ncol = nlevels(sectors) ) # coefficients
e = matrix(0, nrow = nrow(R), ncol = ncol(R) ) # residuals
for( i in 1:nrow(R)){
tmp = lm( R[i,] ~ b - 1 )
f[i,] = tmp$coefficients
e[i,] = tmp$residuals
}
(a) [5 points]
Create a scatter-plot of the values of the two industry factors. Do the factors look independent? Explain your answer based on what you know about asset behaviour.
plot( f[,1], f[,2] )
The two factors are positively correlated. This is expected, since we know that all asset returns have positive correlation, through their dependence on the overall “market” behavior.
(b) [3 points]
Use the sample covariance of the returns ( \(\mathrm{R}\) ) to calculate the minimum-variance portfolio weights.
S = var(R)
Si = solve(S)
w_mv_sample = rowSums(Si)/sum(Si)
(c) [10 points]
Use the industry factor model to estimate the return covariance matrix (from \(\boldsymbol{\beta} \boldsymbol{\Sigma}_F \boldsymbol{\beta}^{\top}+\boldsymbol{S}_e\), and use that to calculate the minimum-variance-portfolio weights.
S_e = diag( apply( e, 2, var ) )
S_f = cov(f)
S_model = S_e + b %*% S_f %*% t(b)
Si = solve( S_model )
w_mv_model = rowSums(Si)/sum(Si)
(d) [3 points]
Use the weights you found in the last two parts to create the returns of each portfolio, and calculate their Sharpe ratio assuming a risk-free rate of 0.
R_mv_sample = R %*% w_mv_sample
(Sharpe_mv_sample = mean( R_mv_sample ) / sd( R_mv_sample ))
## [1] 0.1061
R_mv_model = R %*% w_mv_model
(Sharpe_mv_model = mean( R_mv_model ) / sd( R_mv_model ))
## [1] 0.1321
(e) [4 points]
Plot the cumulative gross returns of the two portfolios on the same axes.
plot( cumprod( 1 + R_mv_model ), type ="l", col = 2 ,
ylab = "gross return", xlab = "day")
lines( cumprod( 1 + R_mv_sample ), col = 1 )
legend( "topleft", col=c(1,2), lwd=2, legend = c("Sample covariance", "Industry factor model covariance"))
2020 Q4 Factor Models
Use the following R code to download daily prices of 10 ETFs, from Jan 1, 2018 to Dec 31, 2019
library(zoo)
library(tseries)
tickers = c("DVEM", "EXT", "HYEM", "LTPZ", "SCHP",
"EDV", "SPMB", "TLT", "GOVT")
S=list()
for(i in 1:length(tickers)){
S[[i]] = get.hist.quote(tickers[i], start='2018-01-01', end='2019-12-31', quote='AdjClose', drop = TRUE)
}
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
## time series starts 2018-01-02
## time series ends 2019-12-30
(a) [5 points]
Calculate the log-returns of the ETFs, and plot the price and return series for the first ETF (DVEM).
plot(S[[1]])
r = lapply(S, FUN=function(x){diff(log(x))}) # log return
R = lapply(r, FUN=function(x){exp(x)-1}) # return
plot(R[[1]])
(b) [5 points]
Use factanal() to fit a 2-factor model to the correlation matrix of the returns. Report the factor loadings and idiosyncratic variances of you model.
Rmat = simplify2array(R)
fmod = factanal( Rmat, factors = 2, lower = 0.005)
b = fmod$loadings
v = fmod$uniquenesses
(c) (10 points]
Simulate 250 daily log-returns using a multivariate Normal distribution with parameters given by the sample means and variances of the ETFs, and correlation matrix given by the previous factor model. Calculate and plot the cumulative net-returns of an equally weighted portfolio over the 10 ETFs.
MU = sapply(R, mean)
SD = sapply(R, sd)
VC = ( b%*%t(b) + diag(v) ) * (SD %*% t(SD))
Rsim = mvtnorm::rmvnorm(250, MU, VC)
R_eqwt = rowMeans( exp( apply(Rsim,2,cumsum) ) - 1 )
plot(R_eqwt, type = "l");
2020 Q5 Monte Carlo Simulation
Consider a European chooser option where the holder gets to decide at time \(T_1\) whether the option becomes a call or a put with fixed strike \(K\) and maturity \(T_2>T_1\).
In other words, the holder “chooses” at time \(T_1\) the form of the option payoff: \(\left(S_{T_2}-K\right)_{+}\)for a call, or \(\left(K-S_{T_2}\right)_{+}\)for a put.
Note that because of put-call parity, i.e. \(C\left(S_{T_1}, T_2-T_1, K\right)-P\left(S_{T_1}, T_2-T_1, K\right)=S_{T_1}-K e^{-r\left(T_2-T_1\right)}\), the holder’s optimal decision at time \(T_1\) is straightforward:
- If \(S_{T_1}-K e^{-r\left(T_2-T_1\right)}>0\), they choose the call (b/c it is more valuable, i.e. \(\left.C\left(S_{T_1}, T_2-T_1, K\right)>P\left(S_{T_1}, T_2-T_1, K\right)\right)\)
- If \(S_{T_1}-K e^{-r\left(T_2-T_1\right)}<0\), they choose the put
Let the current price of the underlying asset be \(S_0=100\), the strike price be \(K=100\), the choosing and expiration times be \(T_1=1\) and \(T_2=2\), and the risk-free rate be \(r=5 \%\), and assume the standard geometric Brownian motion (GBM) asset price dynamics: \[ d S_t=r S_t d t+\sigma S_t d W_t \] with volatility \(\sigma=20 \%\).
(a) [5 points]
Use fExoticOptions: : SimpleChooserOption () to find the exact price of the option.
# library(fExoticOptions)
# SimpleChooserOption(S = S0, X = K, time1 = 1, Time2 = 2,
# r = r, b = r, sigma = v)@price
(b) [10 points]
Perform a simulation with \(n=10,000\) paths for pricing the chooser option. For each path, generate two prices: \(S_{T_1}\) at time \(T_1\) for determining the form of the payoff, and \(S_{T_2}\) at time \(T_2\) for determining the value of the payoff. Report the estimated price and its standard deviation.
set.seed(123);
S0 = K = 100; T1 = 1; T2 = 2; r = .05; v = .2
n = 10000; dT = T2-T1
Z1 = rnorm(10000); Z2 = rnorm(10000)
S1 = S0*exp( (r-v^2/2)*T1 + v*sqrt(T1)*Z1 )
S2 = S1*exp( (r-v^2/2)*dT + v*sqrt(dT)*Z2 )
payoff.b = exp(-r*T2) * ifelse( S1 > K * exp(-r*dT),
pmax( S2 - K, 0),
pmax( K - S2, 0) )
mean(payoff.b)
## [1] 19.97
sd(payoff.b)/sqrt(n)
## [1] 0.2165
(c) [5 points]
Repeat the simulation experiment, but this time use the the Black Scholes formula to value the chosen option at time \(T_1\). In other words, you don’t have to simulate \(S_{T_2}\); just simulate \(S_{T_1}\) and calculate the exact value of the chosen call or put option at time \(T_1\) using Black-Scholes. Report the estimated price and its standard deviation (should be smaller). (Hint: use f0ptions : : GBSOption () to get the Black-Scholes price.)
# library(fOptions)
# payoff.c = exp(-r*T1) * ifelse( S1 > K * exp(-r*dT),
# GBSOption( "c", S1, K, dT, r, r, v )@price,
# GBSOption( "p", S1, K, dT, r, r, v )@price )
# mean(payoff.c)
# sd(payoff.c)/sqrt(N)
(d) [10 points]
Repeat part (b), using the following dynamics for the asset: \[ d S_t=r S_t d t+\sigma \log \left(S_t\right) d W_t \] Note that the price process no longer follows GBM, so use path discretization with \(m=20\) steps per path to approximate \(S_{T_1}\) and \(S_{T_2}\). Report the estimated price and its standard deviation. (Note: the option choice rule at \(T_1\) does not change with asset price dynamics, since put-call parity holds for any model.)
m = 20; Dt = T2/m
Z = matrix( rnorm( n * m ), n, m)
S = matrix( S0, n, m + 1 )
# iterate over all paths
for(i in 1:m){
S[,i+1] = S[,i] + r * S[,i] * Dt + v * log(S[,i]) * sqrt(Dt) * Z[,i]
}
payoff.d = exp(-r*T2) * ifelse( S[, m/2 + 1] > K * exp(-r*dT),
pmax( S[, m + 1] - K, 0),
pmax( K - S[, m + 1], 0) )
mean(payoff.d)
## [1] 9.519
sd(payoff.d) / sqrt(n)
## [1] 0.01262